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

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

         (*.. 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 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);

   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

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

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

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

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

      yields TODO  

   *)

   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

   fun majority_of_coefficients_is_negative 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 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 =

     let

       val (n, q) = fact_quot (lcoeff_up p1') (lcoeff_up p2);

       val remd = drop_tl_zeros ((p1' %* n) %-% (mult_to_deg (deg_up p1') p2 %* q));

       val zeros = for_quot_regarding p1 p2 p1' quot remd

       val _ = writeln_trace p1 p1' p2 quot n q zeros remd ns

       val quot = (replicate zeros 0) @ [q] @ (quot %* n)

       (*div by hand uses q; * n multiplies the previous quot*)

       val ns = n * ns

     in

       if deg_up remd >= deg_up p2

       then div_int_up p1 remd p2 quot ns

       else (ns, quot, remd)

     end

   fun p1 %*/% p2 = 

     let 

       val (n, q, r) = div_int_up p1 p1 p2 [] 1

       val _ = if (! trace_div) then

         (writeln  "  div   .....................................................................";

          writeln ("  div   " ^ PolyML.makestring n ^ " * " ^ PolyML.makestring p1 ^ " = " ^ 

          PolyML.makestring q ^ " ** " ^ PolyML.makestring p2 ^ " ++ " ^ PolyML.makestring r);

          writeln  "  div   ....................................................................."

         ) else ()

     in (n, q, r) end;

 *}

 text {* Here is the check of the algorithm with the values from above: *}

 ML {*

   trace_div := true;

   trace_div_invariant := false;

   val p1 = [2, 2, 2, 2, 2];

   val p2 = [1, 2, 3];

   val (n (* = 27*), 

        q (* = [8, 6, 18]*), 

        r (* = [46, 32]*)) = p1 %*/% p2;

   if (p1 %* n) = ((q %*% p2) %+% r) then () else error "fun %*/% changed"

 *} 

 text {* However, intermediate values increase exponentially: *}

 ML {* (* example from [1].p.84 *)

   val p1 = [~5,2,8,~3,~3,0,1,0,1];

   val p2 = [21,~9,~4,5,0,3];

   val (n (* = 27*), 

        q (* = [22, ~6, ~6, 9]*), 

        r (* = [~597, 378, 376, ~458, 6]*)) = p1 %*/% p2;

   if (p1 %* n) = ((q %*% p2) %+% r) then () else error "fun %*/% changed 1"

 *} 

 ML {* (* example used to develop algorithm by hand *)

   val p1 = [2,2,2, 2,2,2];

   val p2 = [7,8,6];

   val (n (* = 162*), 

        q (* = [55, 15, ~18, 54]*), 

        r (* = [~61, ~221]*)) = p1 %*/% p2;

   if (p1 %* n) = ((q %*% p2) %+% r) then () else error "fun %*/% changed 2"

 *} 

 subsubsection {* The Euclidean Algorithm really works with polynomials? *}

 text {* 

  The generalized division implemented above really can replace the

   division of integers in the Eucliden Algorithm:

 *}

 ML {*

   val trace_Euclid = Unsynchronized.ref true;

   fun writeln_trace_Euclid (a: int list) (b: int list) (n: int) (q: int list) (r: int list) =

     if (! trace_Euclid) then

       writeln ("Euclid "  ^ PolyML.makestring a ^ "  *  "  ^ PolyML.makestring n ^ 

         "  ==  "  ^ PolyML.makestring q  ^ 

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

     else ();

   fun EUCLID_naive_up a [] = primitive_up a

     | EUCLID_naive_up a b =

       if deg_up a < deg_up b then EUCLID_naive_up b a else

       let

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

         val _ = writeln_trace_Euclid a b n q r

       in EUCLID_naive_up b r : unipoly end;

 *}

 ML {*

   (*The example used in the inital mail to Tobias Nipkow made univariate: y -> 1*)

   trace_div := false;

   trace_div_invariant := false;

   trace_Euclid := true;

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

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

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

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

   EUCLID_naive_up a b; (* = [1, ~1]*)          (*( 1 - x) *)

 *}

 ML {*

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

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

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

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

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

 *} 

 ML {* (*NOW this works with generalized division*)

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

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

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

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

   EUCLID_naive_up a b; (* = [~1, 1]*)

 *} 

 ML {* 

   (* from [1].p.84; the generated coefficients are different for an unknown reason *)

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

   val a = [~5,2,8,~3,~3,0,1,0,1]: unipoly;

   val b = [21,~9,~4,5,0,3]: unipoly;

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

   EUCLID_naive_up a b; (* = [1]*)

 *} 

 ML {*

   (* from [1].p.94 *)

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

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

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

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

   EUCLID_naive_up a b; (* = [~2, ~1, 1]*)

 *}

 text {*

   Note: The above algorithm still needs to take into account the case, where 

   constants can be factored out from polynomials, see [1].p.82-84.  

 *} 

 subsubsection {* Interactivity on Euclids Algorithm *}

 text {*

   Euclids Algorithm in modern algebraic notation explains itself if presented

   this way:

 *}

 ML {*

   trace_div := false;

   trace_div_invariant := false;

   trace_Euclid := true;

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

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

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

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

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

 (*

   Euclid [~1, 0, 0, 1]  *  1  ==  [1, 1]  **  [0, ~1, 1]  ++  [~1, 1] 

   Euclid [0, ~1, 1]  *  1  ==  [0, 1]  **  [~1, 1]  ++  [] 

 *)

 *}

 ML {*

   trace_Euclid := true;

   (* from [1].p.94 *)

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

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

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

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

   EUCLID_naive_up a b; (* = [~2, ~1, 1]*)

 (*

 Euclid [~18, ~15, ~20, 12, 20, ~13, 2]  *  1  ==  [1]  **  [8, 28, 22, ~11, ~14, 1, 2]  ++  [~26, ~43, ~42, 23, 34, ~14] 

 Euclid [8, 28, 22, ~11, ~14, 1, 2]  *  98  ==  [~41, ~14]  **  [~26, ~43, ~42, 23, 34, ~14]  ++  [~282, 617, ~168, ~723, 344] 

 Euclid [~26, ~43, ~42, 23, 34, ~14]  *  59168  ==  [787, ~2408]  **  [~282, 617, ~168, ~723, 344]  ++  [~1316434, ~3708859, ~867104, 1525321] 

 Euclid [~282, 617, ~168, ~723, 344]  *  6783102487  ==  [~2345549, 1529768]  **  [~1316434, ~3708859, ~867104, 1525321]  ++  [~5000595353600, ~2500297676800, 2500297676800] 

 Euclid [~1316434, ~3708859, ~867104, 1525321]  *  7289497600  ==  [1919, 4447]  **  [~5000595353600, ~2500297676800, 2500297676800]  ++  [] 

 *)

 *}

 text {*

   A student only needs to observe, how the second factor on the right-hand side

   shifts to the left-hand side; and he or she has to observe the remainder []

   in the last line and take the preceeding remainder as the gcd.

   In order to see that this recursion really computes the gcd, one needs to

   know the fixpoint (or loop invariant):

 *}

 (*

 lemma "a = (q *\<^isub>P b +\<^isub>P r) *\<^isub>s n  \<and>  degree b > degree r  \<longrightarrow>  gcd a b = gcd b r"

 sorry

 *)

 text {*

   The above lemma is analogous to 

   @{thm gcd_add_mult_int} = "gcd ?m (?k * ?m + ?n) = gcd ?m ?n" from HOL/GCD.thy

 *}

 subsubsection {* Interactivity and/or understanding of polynomial division *}

 text {*

   Division raises two requirements wrt. interactivity: 

   (1) a student wants to roughly understand what happens in the algorithm

   (2) a student wants to reproduce the algorithm by hand.

   Case (1) again calls for the fixpoint:

 *}

 ML {*

   trace_div := false;

   trace_div_invariant := true;

   val p1 = [2, 2, 2, 2, 2];

   val p2 = [1, 2, 3];

   val (n (* = 27*), 

        q (* = [8, 6, 18]*), 

        r (* = [46, 32]*)) = p1 %*/% p2;

 (*

   div   [2, 2, 2, 2, 2] * 3 = [0, 0, 2] ** [1, 2, 3] ++ [6, 6, 4, 2, 0] 

   div   [2, 2, 2, 2, 2] * 9 = [0, 2, 6] ** [1, 2, 3] ++ [18, 16, 8, 0, 0] 

   div   [2, 2, 2, 2, 2] * 27 = [8, 6, 18] ** [1, 2, 3] ++ [46, 32, 0, 0, 0] 

 *)

 *}

 text {*

   The trace shows how the quotient, i.e. the first factor on the right-hand side,

   is being increased step by step until the remainder has a degree lower than the

   divisor.

   This representation, however, is inappropriate for reproduction by hand.

   Dividing polynomials by hand is taught in this format:

 *}

 ML {*

   trace_div := true;

   trace_div_invariant := false;

   val p1 = [2, 2, 2, 2, 2];

   val p2 = [1, 2, 3];

   val (n (* = 27*), 

        q (* = [8, 6, 18]*), 

        r (* = [46, 32]*)) = p1 %*/% p2;

 (*

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

   div   ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

   div   [6, 6, 6, 6, 6] 

   div-- [0, 0, 2, 4, 6] <--- [1, 2, 3] * 2 

   div   --------------------------- 

   div   [6, 6, 4, 2, 0]                            quot = [] @ [2] @ (3 * []) 

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

   div   ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

   div   [18, 18, 12, 6] 

   div-- [0, 2, 4, 6] <--- [1, 2, 3] * 2 

   div   --------------------------- 

   div   [18, 16, 8, 0]                            quot = [] @ [2] @ (3 * [2]) 

   div   [18, 16, 8] * 3 

   div   ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

   div   [54, 48, 24] 

   div-- [8, 16, 24] <--- [1, 2, 3] * 8 

   div   --------------------------- 

   div   [46, 32, 0]                            quot = [] @ [8] @ (3 * [2, 6]) 

   div   ..................................................................... 

   div   27 * [2, 2, 2, 2, 2] = [8, 6, 18] ** [1, 2, 3] ++ [46, 32] 

   div   ..................................................................... 

 *)

 *}

 text {*

   The above algorithm is unusual in that the dividend first is multiplied

   such that the quotient can have integer coefficients ("generalised division"). 

   Thus the quotient needs to be multiplied subsequently as well. 

   The algorithm also has to handle cases, where zeros occur in the quotient: 

 *}                      

 ML {*

   trace_div := true;

   val p1 = [3,2,2,2,1,1,1,1];

   val p2 = [1,1,1,1];

   val (n (* = 1*), 

        q (* = [2, 0, 0, 0, 1]*), 

        r (* = [1]*)) = p1 %*/% p2;

 (*

   div   [3, 2, 2, 2, 1, 1, 1, 1] * 1 

   div   ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

   div   [3, 2, 2, 2, 1, 1, 1, 1] 

   div-- [0, 0, 0, 0, 1, 1, 1, 1] <--- [1, 1, 1, 1] * 1 

   div   --------------------------- 

   div   [3, 2, 2, 2, 0, 0, 0, 0]                            quot = [0, 0, 0] @ [1] @ (1 * []) 

   div   [3, 2, 2, 2] * 1 

   div   ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

   div   [3, 2, 2, 2] 

   div-- [2, 2, 2, 2] <--- [1, 1, 1, 1] * 2 

   div   --------------------------- 

   div   [1, 0, 0, 0]                            quot = [] @ [2] @ (1 * [0, 0, 0, 1]) 

   div   ..................................................................... 

   div   1 * [3, 2, 2, 2, 1, 1, 1, 1] = [2, 0, 0, 0, 1] ** [1, 1, 1, 1] ++ [1] 

   div   ..................................................................... 

 *)

 *}

 subsubsection {*  *}

 text {*

 *}

 ML {*

 *}

 subsection {* Technical requirements raised by pedagogical questions *}

 subsubsection {* Calculations might be different from interpretation *}

 text {*

   Until now the demonstration cases for Lucas-Interpretation were

   simplification, equation solving and ad-hoc algorithms in Signal Processing

   and Structural Engineering. In these cases execution of the functional program 

   coincided with the calculation generated as a side-effect of interpretation:

   The tactics "Take", "Rewrite", "Rewrite_Set" and "Substitute", when executed,

   immediately created respective output and handed over control to the user.

   Appropriate user-input finally created a "calculation" very close to 

   what would be written on paper.

   Now, in Computer Algebra, the programs contain lots of sophisticated functions.

   Simple traces, for instance arguments and values of function calls, comprise too 

   many data to be intuitively comprehensible.

   The examples in the case study suggest to have mechanisms for creating

   steps for calculations, which are somehow abstracted from pure execution.

   Interestingly these abstractions are close to the fixpoints/invariants

   of the respective recursive functions --- see 

   section "Interactivity on Euclids Algorithm" and section "Interactivity and/or 

   understanding of polynomial division" above.

   Such mechanisms might be specific tactics, which take data from the 

   local environment and the local context for presentation to the student.

   Thus one such tactic is related to several lines of code in a function.

   These tactics would not contribute to automated execution as done by "value",

   for instance.

   However, the effect of these tactics for resuming execution upon completion 

   of input needs to be clarified.

 *}

 subsubsection {* One function might generate various representations *}

 text {*

   The section "Interactivity and/or understanding of polynomial division" suggests

   to have at least two kinds of representation for "calculations": 

   (1) a student wants to roughly understand what happens in the algorithm

   (2) a student wants to reproduce the algorithm by hand.

   A solution to this requirement of variety is to create related functions 

   with different combinations of the specific tactics as mentioned above.

   These functions can be selected by the dialogue module. The computational

   effect of such related functions is the same, the interactive behaviour

   is different (much more different as variations of "Rewrite", for instance.)

 *}

 subsubsection {* One "specific tactic" might create several lines of calculation *}

 text {*

   The case of "division of polynomial" suggests to have one "specific tactic"

   as mentioned above, which has one location in the function (e.g. "writeln_trace"

   in "div_int_up"), but creates several lines of a calculation; see the example in

   section "Can division of polynomials be adapted":

         [  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

       :

       :

   Interaction on these lines should resemble calculation with paper and pencil.

 *}

 subsubsection {*  *}

 text {*

 *}

 ML {*

 *}

 end