header {* GCD for polynomials, implemented using the function package (_FP) *}

 theory GCD_Poly_FP 

 imports "~~/src/HOL/Algebra/poly/Polynomial" (*imports ~~/src/HOL/Library/Polynomial ...2012*)

         "~~/src/HOL/Number_Theory/Primes"

 (* WN130304.TODO see Algebra/poly/LongDiv.thy: lcoeff_def*)

 begin

 text {* 

   This code has been translated from GCD_Poly.thy by Diana Meindl,

   who follows Franz Winkler, Polynomial algorithyms in computer algebra, Springer 1996.

   Winkler's original identifiers are in test/./gcd_poly_winkler.sml;

   test/../gcd_poly.sml documents the changes towards GCD_Poly.thy;

   the style of GCD_Poly.thy has been adapted to the function package.

 *}

 section {* gcd for univariate polynomials *}

 type_synonym unipoly = "int list" (*TODO: compare Polynomial.thy*)

 value "[0, 1, 2, 3, 4, 5] :: unipoly"

 subsection {* auxiliary functions *}

 (* a variant for div: 

   5 div 2 = 2; ~5 div 2 = ~3; BUT WEE NEED ~5 div2 2 = ~2; *)

 definition div2 :: "int \<Rightarrow> int \<Rightarrow> int" (infixl "div2" 70) where

 "a div2 b = (if a div b < 0 then (\<bar>a\<bar> div \<bar>b\<bar>) * -1 else a div b)"

 value " 5 div2  2 =  2"

 value "-5 div2  2 = -2"

 value "-5 div2 -2 =  2"

 value " 5 div2 -2 = -2"

 value "gcd (15::int) (6::int) = 3"

 value "gcd (10::int) (3::int) = 1"

 value "gcd (6::int) (24::int) = 6"

 (* drop tail elements equal 0 *)

 primrec drop_hd_zeros :: "int list \<Rightarrow> int list" where

 "drop_hd_zeros (p # ps) = (if p = 0 then drop_hd_zeros ps else (p # ps))"

 definition drop_tl_zeros :: "int list \<Rightarrow> int list" where

 "drop_tl_zeros p = (rev o drop_hd_zeros o rev) p"

 subsection {* modulo calculations for integers *}

 (* modi is just local for mod_inv *)

 function modi :: "int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat" where

 "modi r rold m rinv = 

 (if m \<le> rinv

   then 0 else

     if r mod (int m) = 1 

     then rinv 

     else modi (rold * ((int rinv) + 1)) rold m (rinv + 1))"

 by auto

 termination by (relation "measure (\<lambda>(r, rold, m, rinv). m - rinv)") auto

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

    mod_inv r m = s

    assumes True

    yields r * s mod m = 1 *)

 definition mod_inv :: "int \<Rightarrow> nat \<Rightarrow> nat" where "mod_inv r m = modi r r m 1"

 value "modi 5 5 7 1   = 3"

 value "modi 3 3 7 1   = 5"

 value "modi 4 4 339 1 = 85"

 value "mod_inv 5 7    = 3"

 value "mod_inv 3 7    = 5"

 value "mod_inv 4 339  = 85"

 value "mod_inv -5 7    = 4"

 value "mod_inv -3 7    = 2"

 value "mod_inv -4 339  = 254"

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

 definition mod_div :: "int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat" where

 "mod_div a b m = (nat a) * (mod_inv b m) mod m"

 (*if mod_div 21 4 5 = 4 then () else error "mod_div 21 4 5 = 4 changed";

 if mod_div 1 4 5 = 4 then () else error "mod_div 1 4 5 = 4 changed";

 if mod_div 0 4 5 = 0 then () else error "mod_div 0 4 5 = 0 changed";

 *)

 value "mod_div 21 3 5 = 2"            value "(21::int) mod 5 = (3 * 2) mod 5"

 (*             21/3   = 7 mod 5               21       mod 5 =    6    mod 5

                       = 2                                  1               1 *)

 value "mod_div 22 3 5 = 4"            value "(22::int) mod 5 = (3 * 4) mod 5"

 (*             22/3   = -------               22       mod 5 =   12    mod 5

                       = 4                                  2               2 *)

 value "mod_div 23 3 5 = 1"            value "(23::int) mod 5 = (3 * 1) mod 5"

 (*             23/3   = -------               23       mod 5 =    3    mod 5

                       = 1                                  3               3 *)

 value "mod_div 24 3 5 = 3"            value "(24::int) mod 5 = (3 * 3) mod 5"

 (*             24/3   = -------               24       mod 5 =    9    mod 5

                       = 3                                  4               4 *)

 value "mod_div 25 3 5 = 0"            value "(25::int) mod 5 = (3 * 0) mod 5"

 (*             25/3   = -------               25       mod 5 =    0    mod 5

                       = 0                                  0               0 *)

 value "mod_div 21 4 5 = 4"            value "(21::int) mod 5 = (4 * 4) mod 5"

 value "mod_div  1 4 5 = 4"            value "( 1::int) mod 5 = (4 * 4) mod 5"

 value "mod_div  0 4 5 = 0"            value "( 0::int) mod 5 = (0 * 4) mod 5"

 (* root1 is just local to approx_root *)

 function root1 :: "int \<Rightarrow> nat \<Rightarrow> nat" where

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

 by auto

 termination by (relation "measure (\<lambda>(a, n). nat (a - (int (n * n))))") auto

 (* required just for one approximation:

    approx_root :: nat \<Rightarrow> nat

    approx_root a = r

    assumes True

    yields r * r \<ge> a *)

 definition approx_root :: "int \<Rightarrow> nat" where "approx_root a = root1 a 1"

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

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

    assumes True

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

 definition chinese_remainder :: "int \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat" where

 "chinese_remainder r1 m1 r2 m2 = 

   ((nat (r1 mod (int m1))) + 

     (nat (((r2 - (r1 mod (int m1))) * (int (mod_inv (int m1) m2))) mod (int m2))) * m1)"

 value "chinese_remainder 17 9 3 4  = 35"

 value "chinese_remainder  7 2 6 11 = 17"

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

 (* 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 :: "nat list \<Rightarrow> nat \<Rightarrow> bool" where

 "is_prime ps n = 

 (if List.length ps > 0

 then 

   if (n mod (List.hd ps)) = 0

   then False

   else is_prime (List.tl ps) n

 else True)"

 declare is_prime.simps [simp del] -- "make_primes, next_prime_not_dvd"

 value "is_prime [2, 3] 2  = False"    --"... precondition!"

 value "is_prime [2, 3] 3  = False"    --"... precondition!"

 value "is_prime [2, 3] 4  = False"

 value "is_prime [2, 3] 5  = True"

 value "is_prime [2, 3, 5] 5  = False" --"... precondition!"

 value "is_prime [2, 3] 6  = False"

 value "is_prime [2, 3] 7  = True"

 value "is_prime [2, 3] 25 = True"     -- "... because 5 not in list"

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

 function make_primes :: "nat list \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat list" where

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

 by pat_completeness auto --"simp del: is_prime.simps <-- declare"

 (*Calls:

   a) (ps, last_p, n) ~> (ps @ [last_p + 2], last_p + 2, n)

   b) (ps, last_p, n) ~> (ps,                last_p + 2, n) 

 Measures:

   1) \<lambda>p. size (snd (snd p))     --"n"

   2) \<lambda>p. size (fst (snd p))     --"last_p"

   3) \<lambda>p. size (snd p)           --"(last_p, n)"

   4) \<lambda>p. list_size size (fst p) --"ps"

   5) \<lambda>p. length (fst p)         --"ps"

   6) size

 Result matrix:

     1  2  3  4  5  6 

 a:  <= ?  <= ?  ?  <=

 b:  <= ?  <= <= <= <= *)

 find_theorems "_ - _ < _ = _"

 find_theorems "?n < ?n + _" 

   thm Rings.linordered_semidom_class.less_add_one --"?a < ?a + 1"

 find_theorems "_ + ?a < _ + ?a = _" 

   thm Groups.ordered_ab_semigroup_add_imp_le_class.add_less_cancel_right

 find_theorems "_ - ?a < _ - ?a = _" 

   thm Nat.less_diff_iff

 lemma termin_1: "n - Suc (length ps) < n - length ps" (*? \<not> True ?*) sorry

 termination make_primes (*by lexicographic_order +PROOF:2 GOLAS / size_change LOOPS*)

 apply (relation "measures [\<lambda>(ps, last_p, n). n - (length ps),

                            \<lambda>(ps, last_p, n). n + 2 - last_p]")

 apply auto

 (*\<lambda>(ps, last_p, n). last_p - (length ps):

 goal (4 subgoals):

  1. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> Suc last_p - length ps < last_p - length ps \<Longrightarrow> Suc last_p - length ps = last_p - length ps

  2. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> Suc last_p - length ps < last_p - length ps \<Longrightarrow> n - last_p < Suc (Suc n) - last_p

  3. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> \<not> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> Suc (Suc last_p) - length ps < last_p - length ps \<Longrightarrow> Suc (Suc last_p) - length ps = last_p - length ps

  4. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> \<not> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> Suc (Suc last_p) - length ps < last_p - length ps \<Longrightarrow> n - last_p < Suc (Suc n) - last_p

 \<lambda>(ps, last_p, n). n - (length ps):

 goal (2 subgoals):

  1. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> \<not> n - Suc (length ps) < n - length ps \<Longrightarrow> n - last_p < Suc (Suc n) - last_p

  2. \<And>ps last_p n. \<not> n \<le> last ps \<Longrightarrow> \<not> is_prime ps (Suc (Suc last_p)) \<Longrightarrow> n - last_p < Suc (Suc n) - last_p*)

 sorry

 declare make_primes.simps [simp del] -- "next_prime_not_dvd"

 value "make_primes [2, 3] 3  3 = [2, 3]"

 value "make_primes [2, 3] 3  4 = [2, 3, 5]"

 value "make_primes [2, 3] 3  5 = [2, 3, 5]"

 value "make_primes [2, 3] 3  6 = [2, 3, 5, 7]"

 value "make_primes [2, 3] 3  7 = [2, 3, 5, 7]"

 value "make_primes [2, 3] 3  8 = [2, 3, 5, 7, 11]"

 value "make_primes [2, 3] 3  9 = [2, 3, 5, 7, 11]"

 value "make_primes [2, 3] 3 10 = [2, 3, 5, 7, 11]"

 value "make_primes [2, 3] 3 11 = [2, 3, 5, 7, 11]"

 value "make_primes [2, 3] 3 12 = [2, 3, 5, 7, 11, 13]"

 value "make_primes [2, 3] 3 13 = [2, 3, 5, 7, 11, 13]"

 value "make_primes [2, 3] 3 14 = [2, 3, 5, 7, 11, 13, 17]"

 value "make_primes [2, 3] 3 15 = [2, 3, 5, 7, 11, 13, 17]"

 value "make_primes [2, 3] 3 16 = [2, 3, 5, 7, 11, 13, 17]"

 value "make_primes [2, 3] 3 17 = [2, 3, 5, 7, 11, 13, 17]"

 value "make_primes [2, 3] 3 18 = [2, 3, 5, 7, 11, 13, 17, 19]"

 value "make_primes [2, 3] 3 19 = [2, 3, 5, 7, 11, 13, 17, 19]"

 value "make_primes [2, 3, 5, 7] 7 4 = [2, 3, 5, 7]"

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

 definition primes_upto :: "nat \<Rightarrow> nat list" where

 "primes_upto n = (if n < 3 then [2] else make_primes [2, 3] 3 n)"

 value "primes_upto  0 = [2]"

 value "primes_upto  1 = [2]"

 value "primes_upto  2 = [2]"

 value "primes_upto  3 = [2, 3]"

 value "primes_upto  4 = [2, 3, 5]"

 value "primes_upto  5 = [2, 3, 5]"

 value "primes_upto  6 = [2, 3, 5, 7]"

 value "primes_upto  7 = [2, 3, 5, 7]"

 value "primes_upto  8 = [2, 3, 5, 7, 11]"

 value "primes_upto  9 = [2, 3, 5, 7, 11]"

 value "primes_upto 10 = [2, 3, 5, 7, 11]"

 value "primes_upto 11 = [2, 3, 5, 7, 11]"

 lemma primes_upto_0: "last (primes_upto n) > 0" (*see above*) sorry

 lemma primes_upto_1: "last (primes_upto (Suc n)) > n" (*see above*) sorry

 lemma primes_upto_2: "last (primes_upto n) >= n" (*see above*) sorry

 lemma termin_next_prime_not_dvd:

   shows "last (primes_upto (Suc n1)) * q - last (primes_upto (Suc n1))

        < last (primes_upto (Suc n1)) * q - n1"

 proof -

   from primes_upto_1 have "n1 < last (primes_upto (Suc n1))" by auto

   from this have                "(int n1) - (int (last (primes_upto (Suc n1)) * q))

     < (int (last (primes_upto (Suc n1)))) - (int (last (primes_upto (Suc n1)) * q))" by auto

   from this show ?thesis sorry

 qed (*?FLORIAN : lifting nat -> int ?*)

 find_theorems "_ - _ < _ - _ = _" thm Nat.less_diff_iff

 find_theorems "_ * _ < _ * _ = _" thm Nat.Suc_mult_less_cancel1

 find_theorems "-_ < -_ = _" thm Groups.ordered_ab_group_add_class.neg_less_iff_less

 (* max's' is analogous to Integer.gcds *)

 definition maxs :: "nat list \<Rightarrow> nat" where "maxs ns = List.fold max ns (List.hd ns)"

 value "maxs [5, 3, 7, 1, 2, 4] = 7"

 (* 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  assumes "0 < q" 

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

 function next_prime_not_dvd :: "nat \<Rightarrow> nat \<Rightarrow> nat" (infixl "next'_prime'_not'_dvd" 70) where

 "n1 next_prime_not_dvd n2 = 

 (let

   ps = primes_upto (n1 + 1);

   nxt = last ps

 in

   if n2 mod nxt \<noteq> 0

   then nxt

   else nxt next_prime_not_dvd n2)"

 by auto --"simp del: is_prime.simps, make_primes.simps, primes_upto.simps <-- declare"

 termination (*next_prime_not_dvd: lexicographic_order  +PROOF:?lifting nat-int? / size_change: Failed to apply initial proof method*)

   apply (relation "measure (\<lambda>(n1, n2). n2 - n1)") 

   apply auto

   apply (rule termin_next_prime_not_dvd) (*?FLORIAN: IN Try_GCD... this is not necessary!?!+*)

 done

 value "1 next_prime_not_dvd 15 = 2"

 value "2 next_prime_not_dvd 15 = 7"

 value "3 next_prime_not_dvd 15 = 7"

 value "4 next_prime_not_dvd 15 = 7"

 value "5 next_prime_not_dvd 15 = 7"

 value "6 next_prime_not_dvd 15 = 7"

 value "7 next_prime_not_dvd 15 =11"

 subsection {* basic notions for univariate polynomials *}

 (* not in List.thy, copy from library.ML *)

 fun nth_drop :: "nat \<Rightarrow> 'a list \<Rightarrow> 'a list" where

 "nth_drop n xs = List.take n xs @ List.drop (n + 1) xs"

 value "nth_drop 0 []              = []"

 value "nth_drop 0 [1, 2, 3::int]  = [2, 3]"

 value "nth_drop 2 [1, 2, 3::int]  = [1, 2]"

 value "nth_drop 77 [1, 2, 3::int] = [1, 2, 3]"

 (* leading coefficient *)

 definition lcoeff_up :: "unipoly \<Rightarrow> int" where "lcoeff_up p = (last o drop_tl_zeros) p"

 value "lcoeff_up [3, 4, 5, 6]    = 6"

 value "lcoeff_up [3, 4, 5, 6, 0] = 6"

 (* drop leading coefficients equal 0 *)

 (* THESE MAKE value "gcd_up [-1, 0 ,1] [0, 1, 1] = [1, 1]" LOOP 

   WHILE SML-VERSION WORKS:  (*?FLORIAN*)

 (**)definition drop_lc0_up :: "unipoly \<Rightarrow> unipoly" where 

   "drop_lc0_up p = drop_tl_zeros p"

 (**)fun drop_lc0_up :: "unipoly \<Rightarrow> unipoly" where 

   "drop_lc0_up p = drop_tl_zeros p"

 *)

 fun drop_lc0_up :: "unipoly \<Rightarrow> unipoly" where 

 "drop_lc0_up [] = []" |

 "drop_lc0_up p = 

   (let l = List.length p - 1

   in

     if p ! l = 0

     then drop_lc0_up (nth_drop l p)

     else p)"

 value "drop_lc0_up [0, 1, 2, 3, 4, 5, 0, 0] = [0, 1, 2, 3, 4, 5]"

 value "drop_lc0_up [0, 1, 2, 3, 4, 5]       = [0, 1, 2, 3, 4, 5]"

 (* degree *)

 (* THE VERSIONS BELOW CREATE THE RESULT value "[-18,.." --> value "[-1,.." = [1] 

   ALTHOUGH THE TESTS BELOW SEEM TO WORK AND THE SAME CODE WORKS IN ML ?FLORIAN*)

 (*function deg_up :: "unipoly \<Rightarrow> nat" where

   "deg_up p = ((op - 1) o length o drop_lc0_up) p"

 by auto termination sorry*)

 (*fun deg_up :: "unipoly \<Rightarrow> nat" where

   "deg_up p = ((op - 1) o length o drop_lc0_up) p"*)

 (*definition deg_up :: "unipoly \<Rightarrow> nat" where

   "deg_up p = ((op - 1) o length o drop_lc0_up) p"*)

 function deg_up :: "unipoly \<Rightarrow> nat" where

 "deg_up p =

 (let len = List.length p - 1

 in

   if p ! len \<noteq> 0 then len else deg_up (nth_drop len p))"

 by (*pat_completeness*) auto 

 lemma min_Suc:

   fixes a::nat

   assumes "0 < a"

   shows "min a (a - Suc 0) < a"

 proof -

   from min_def have 1: "min a (a - Suc 0) = a - Suc 0" by auto

   from `0 < a` have 2: "a - Suc 0 < a" by auto

   from 1 2 show ?thesis by auto

 qed

 find_theorems "min _ _ = (_::nat)"

 find_theorems "min _ _ "

   thm Orderings.ord_class.min_def --"min ?a ?b = (if ?a \<le> ?b then ?a else ?b)"

 termination deg_up (*by lexicographic_order  +PROOF:STRANGE GOAL+definition / by size_change: Failed to apply initial proof method*)

   apply (relation "measure (\<lambda>(p). length p)")

   apply auto

 (*..RELATED?: https://lists.cam.ac.uk/mailman/htdig/cl-isabelle-users/2013-May/msg00075.html*)

   apply (rule min_Suc) (* 1. \<And>p. p ! (length p - Suc 0) = 0 \<Longrightarrow> 0 < length p*)

   apply auto     (*?????: 1. [] ! 0 = 0 \<Longrightarrow> False *)

 sorry 

 find_theorems "nth _ _"

   thm List.nth_mem        (*?n < length ?xs \<Longrightarrow> ?xs ! ?n \<in> set ?xs*)

           (*(length p - Suc 0) < length ?xs *)

   thm   List.set_conv_nth (*set ?xs = {?xs ! i |i. i < length ?xs}*)

 find_theorems "set _"    (*found 378 theorem(s)*)

 find_theorems "length _" (*found 241 theorem(s)*)

 value "[1,2,3,4,5::int] ! 2"

 value "[1,2,3,4,5::int] ! 4"

 value "[1::int] ! 0"

 value "([]::int list) ! 0"

 (*Calls:

   a) p ~> nth_drop x p

 Measures:

   1) list_size (nat \<circ> abs)

   2) length

 Result matrix:

     1  2 

 a:  ?  <= *)

 value "deg_up [3, 4, 5, 6]    = 3"

 value "deg_up [3, 4, 5, 6, 0] = 3"

 value "deg_up [1, 0, 3, 0, 0] = 2"

 (* norm is just local to normalise *)

 fun norm :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> unipoly" where

 "norm p nrm m lcp i = 

 (if i = 0

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

   else norm p ([int (mod_div (p ! i) lcp m)] @ nrm) m lcp (i - 1))"

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

      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 :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where

 "normalise p m = 

 (let

  p = drop_lc0_up p;

  lcp = lcoeff_up p

 in 

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

 declare normalise.simps [simp del] --"HENSEL_lifting_up"

 value "normalise [-18, -15, -20, 12, 20, -13, 2] 5 = [1, 0, 0, 1, 0, 1, 1]"

 value "normalise [9, 6, 3] 10                      = [3, 2, 1]"

 subsection {* addition, multiplication, division *}

 (* scalar multiplication *)

 definition  mult_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "%*" 70) where

 "p %* m = List.map (op * m) p"

 value "[5, 4, 7, 8, 1] %* 5   = [25, 20, 35, 40, 5]"

 value "[5, 4, -7, 8, -1] %* 5 = [25, 20, -35, 40, -5]"

 (* scalar divison *)

 definition swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where "swapf f a b = f b a"

 definition div_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "div'_ups" (* %/ error FLORIAN*) 70) where

 "p div_ups m = map (swapf op div2 m) p"

 value "[4, 3, 2, 5, 6] div_ups 3 = [1, 1, 0, 1, 2]"

 value "[4, 3, 2, 0] div_ups 3    = [1, 1, 0, 0]"

 (* not in List.thy, copy from library.ML *)

 fun map2 :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'a list \<Rightarrow> 'b list \<Rightarrow> 'c list" where

 "map2 _ [] [] = []" |

 "map2 f (x # xs) (y # ys) = f x y # map2 f xs ys" |

 "map2 _ _ _ = []" (*raise ListPair.UnequalLengths*)

 (* minus of polys *)

 definition minus_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" (infixl "%-%" 70) where

 "p1 %-% p2 = map2 (op -) p1 p2"

 value "[8, -7, 0, 1] %-% [-2, 2, 3, 0]     = [10, -9, -3, 1]"

 value "[8, 7, 6, 5, 4] %-% [2, 2, 3, 1, 1] = [6, 5, 3, 4, 3]"

 (* lemmata for pattern compatibility in dvd_up *)

 lemma ex_falso_1: "([d], [p]) = (v # vb # vc, ps) \<Longrightarrow> QUODLIBET" by simp

 lemma ex_falso_2: "([], ps) = (v # vb # vc, psa) \<Longrightarrow> QUODLIBET" by simp

 lemma ex_falso_3: "(ds, []) = (dsa, v # vb # vc) \<Longrightarrow> QUODLIBET" by simp

 function (sequential) dvd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> bool" (infixl "%|%" 70) where

 "[d] %|% [p] = ((\<bar>d\<bar> \<le> \<bar>p\<bar>) \<and> (p mod d = 0))" |

 "ds %|% ps =

   (let 

     ds = drop_lc0_up ds; ps = drop_lc0_up ps;

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

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

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

   in 

     if rest = [] then True else 

       if quot \<noteq> 0 \<and> List.length ds \<le> List.length rest then ds %|% rest else False)"

 apply pat_completeness

 apply simp

 apply simp

 defer (*a "mixed" obligation*)

 apply simp

 defer (*a "mixed" obligation*)

 apply simp

 defer (*a "mixed" obligation*)

 apply simp

 apply simp

 apply simp (* > 1 sec IMPROVED BY declare simp del: 

   centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)

 apply simp

 apply simp (* > 1 sec IMPROVED BY declare simp del: 

   centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)

 apply simp

 defer (*a "mixed" obligation*)

 apply simp (* > 1 sec IMPROVED BY declare simp del: 

   centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)

 defer (*a "mixed" obligation*)

 apply (rule ex_falso_1)

 apply simp

 apply (rule ex_falso_2)

 apply simp

 apply (rule ex_falso_3)

 apply simp

 apply (rule ex_falso_1)

 apply simp

 done (* > 1 sec IMPROVED BY declare simp del: 

   centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)

 termination (*dvd_up: by lexicographic_order LOOPS  +PROOF:5 HUGE GOALS/ size_change LOOPS*)

 using [[linarith_split_limit = 999]]

 apply (relation "measure (\<lambda>(ds, ps). length (let 

                       ds = drop_lc0_up ds; ps = drop_lc0_up ps;

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

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

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

                     in 

                       ds) - length ds)")

 apply auto

 (*apply (relation "measure (\<lambda>(ds, ps). length ps - length ds)"):

  1. wf (measure (\<lambda>(ds, ps). length ps - length ds))

  2. \<And>ps x xa xb xc xd. x = drop_lc0_up [] \<Longrightarrow> xa = drop_lc0_up ps \<Longrightarrow> xb = replicate (length xa - length x) 0 @ x \<Longrightarrow> xc = lcoeff_up xa div2 lcoeff_up xb \<Longrightarrow> xd = drop_lc0_up (xa %-% (xb %* xc)) \<Longrightarrow> xd \<noteq> [] \<Longrightarrow> xc \<noteq> 0 \<and> length x \<le> length xd \<Longrightarrow> ((x, xd), [], ps) \<in> measure (\<lambda>(ds, ps). length ps - length ds)

  3. \<And>v vb vc ps x xa xb xc xd. x = drop_lc0_up (v # vb # vc) \<Longrightarrow> xa = drop_lc0_up ps \<Longrightarrow> xb = replicate (length xa - length x) 0 @ x \<Longrightarrow> xc = lcoeff_up xa div2 lcoeff_up xb \<Longrightarrow> xd = drop_lc0_up (xa %-% (xb %* xc)) \<Longrightarrow> xd \<noteq> [] \<Longrightarrow> xc \<noteq> 0 \<and> length x \<le> length xd \<Longrightarrow> ((x, xd), v # vb # vc, ps) \<in> measure (\<lambda>(ds, ps). length ps - length ds)

  4. \<And>ds x xa xb xc xd. x = drop_lc0_up ds \<Longrightarrow> xa = drop_lc0_up [] \<Longrightarrow> xb = replicate (length xa - length x) 0 @ x \<Longrightarrow> xc = lcoeff_up xa div2 lcoeff_up xb \<Longrightarrow> xd = drop_lc0_up (xa %-% (xb %* xc)) \<Longrightarrow> xd \<noteq> [] \<Longrightarrow> xc \<noteq> 0 \<and> length x \<le> length xd \<Longrightarrow> ((x, xd), ds, []) \<in> measure (\<lambda>(ds, ps). length ps - length ds)

  5. \<And>ds v vb vc x xa xb xc xd. x = drop_lc0_up ds \<Longrightarrow> xa = drop_lc0_up (v # vb # vc) \<Longrightarrow> xb = replicate (length xa - length x) 0 @ x \<Longrightarrow> xc = lcoeff_up xa div2 lcoeff_up xb \<Longrightarrow> xd = drop_lc0_up (xa %-% (xb %* xc)) \<Longrightarrow> xd \<noteq> [] \<Longrightarrow> xc \<noteq> 0 \<and> length x \<le> length xd \<Longrightarrow> ((x, xd), ds, v # vb # vc) \<in> measure (\<lambda>(ds, ps). length ps - length ds)

 apply auto

  6 subgoals, even longer

 *)

 sorry 

 value "[4] %|% [6]                 = False"

 value "[8] %|% [16, 0]             = True"

 value "[3, 2] %|% [0, 0, 9, 12, 4] = True"

 value "[8, 0] %|% [16]             = True"

 subsection {* normalisation and Landau-Mignotte bound *}

 (* centr is just local to centr_up *)

 definition centr :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int" where

 "centr m mid p_i = (if mid < p_i then p_i - (int m) else p_i)"

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

 definition centr_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where

 "centr_up p m =

 (let

   mi = (int m) div2 2;

   mid = if (int m) mod mi = 0 then mi else mi + 1

 in map (centr m mid) p)"

 value "centr_up [7, 3, 5, 8, 1, 3] 10 = [-3, 3, 5, -2, 1, 3]"

 value "centr_up [1, 2, 3, 4, 5] 2     = [1, 0, 1, 2, 3]"

 value "centr_up [1, 2, 3, 4, 5] 3     = [1, -1, 0, 1, 2]"

 value "centr_up [1, 2, 3, 4, 5] 4     = [1, 2, -1, 0, 1]"

 value "centr_up [1, 2, 3, 4, 5] 5     = [1, 2, 3, -1, 0]"

 value "centr_up [1, 2, 3, 4, 5] 6     = [1, 2, 3, -2, -1]"

 value "centr_up [1, 2, 3, 4, 5] 7     = [1, 2, 3, 4, -2]"

 value "centr_up [1, 2, 3, 4, 5] 8     = [1, 2, 3, 4, -3]"

 value "centr_up [1, 2, 3, 4, 5] 9     = [1, 2, 3, 4, 5]"

 value "centr_up [1, 2, 3, 4, 5] 10    = [1, 2, 3, 4, 5]"

 value "centr_up [1, 2, 3, 4, 5] 11    = [1, 2, 3, 4, 5]"

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

 definition sum_lmb :: "unipoly \<Rightarrow> nat \<Rightarrow> int" where

 "sum_lmb p exp = List.fold ((op +) o (swapf power exp)) p 0"

 value "sum_lmb [-1, 2, -3, 4, -5] 1 = -3"

 value "sum_lmb [-1, 2, -3, 4, -5] 2 = 55"

 value "sum_lmb [-1, 2, -3, 4, -5] 3 = -81"

 value "sum_lmb [-1, 2, -3, 4, -5] 4 = 979"

 value "sum_lmb [-1, 2, -3, 4, -5] 5 = -2313"

 value "sum_lmb [-1, 2, -3, 4, -5] 6 = 20515"

 (* LANDAU_MIGNOTTE_bound :: centr_up \<Rightarrow> unipoly => 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)*)

 definition LANDAU_MIGNOTTE_bound :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat" where

 "LANDAU_MIGNOTTE_bound p1 p2 =

   ((power 2 (min (deg_up p1) (deg_up p2))) * (nat (gcd (lcoeff_up p1) (lcoeff_up p2))) * 

   (nat (min (abs ((int (approx_root (sum_lmb p1 2))) div2 -(lcoeff_up p1)))

             (abs ((int (approx_root (sum_lmb p2 2))) div2 -(lcoeff_up p2))))))"

 value "LANDAU_MIGNOTTE_bound [1] [4, 5]          = 1"

 value "LANDAU_MIGNOTTE_bound [1, 2] [4, 5]       = 2"

 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5]    = 2"

 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4]       = 1"

 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5]    = 2"

 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5, 6] = 12"

 value "LANDAU_MIGNOTTE_bound [-1] [4, 5]            = 1"

 value "LANDAU_MIGNOTTE_bound [-1, 2] [4, 5]         = 2"

 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5]    = 2"

 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4]        = 1"

 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5]    = 2"

 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5, 6] = 12"

 subsection {* modulo calculations for polynomials *}

 (* pair is just local to chinese_remainder_up, is "op ~~" in library.ML *)

 fun pair :: "unipoly \<Rightarrow> unipoly \<Rightarrow> ((int \<times> int) list)" (infix "pair" 4) where

 "([] pair []) = []" |

 "([] pair ys) = []" | (*raise ListPair.UnequalLengths*)

 "(xs pair []) = []" | (*raise ListPair.UnequalLengths*)

 "((x#xs) pair (y#ys)) = (x, y) # (xs pair ys)"

 fun chinese_rem :: "nat \<times> nat \<Rightarrow> int \<times> int \<Rightarrow> int" where

 "chinese_rem (m1, m2) (p1, p2) = (int (chinese_remainder p1 m1 p2 m2))"

   (* 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 :: "nat \<times> nat \<Rightarrow> unipoly \<times> unipoly \<Rightarrow> unipoly" where

 "chinese_remainder_up (m1, m2) (p1, p2) = map (chinese_rem (m1, m2)) (p1 pair p2)"

 value "chinese_remainder_up (5, 7) ([2, 2, 4, 3], [3, 2, 3, 5]) = [17, 2, 24, 33]"

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

 definition mod' :: "nat \<Rightarrow> int \<Rightarrow> int" where "mod' m i = i mod (int m)"

 definition mod_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" (infixl "mod'_up" 70) where

 "p mod_up m = map (mod' m) p"

 value "[5, 4, 7, 8, 1] mod_up 5 = [0, 4, 2, 3, 1]"

 value "[5, 4,-7, 8,-1] mod_up 5 = [0, 4, 3, 3, 4]"

 (* euclidean algoritm in Z_p[x/m].

    mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly

    mod_up_gcd p1 p2 m = g

      assumes 

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

 function mod_up_gcd :: "unipoly \<Rightarrow> unipoly  \<Rightarrow> nat \<Rightarrow> unipoly" where

 "mod_up_gcd p1 p2 m = 

 (let 

   p1m = p1 mod_up m;

   p2m = drop_lc0_up (p2 mod_up m);

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

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

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

 in 

   if rest = [] then p2 else

     if List.length rest < List.length p2

     then mod_up_gcd p2 rest m 

     else mod_up_gcd rest p2 m)"

 by auto

 termination mod_up_gcd (*by lexicographic_order +PROOF:3 HUGE SUGOALS / size_change LOOPS*)

 apply (relation "measures [\<lambda>(p1, p2, m). length p2 - length (let 

                                 p1m = p1 mod_up m;

                                 p2m = drop_lc0_up (p2 mod_up m);

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

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

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

                               in rest),

                            \<lambda>(p1, p2, m). length p1 - length (let 

                                 p1m = p1 mod_up m;

                                 p2m = drop_lc0_up (p2 mod_up m);

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

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

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

                               in rest)]")

 (*apply auto ..LOOPS*)

 sorry

 (*Calls:

   a) (p1, p2, m) ~> (p2, xab, m)

   b) (p1, p2, m) ~> (xab, p2, m)

 Measures:

   1) \<lambda>p. size (snd (snd p))                   m

   2) \<lambda>p. list_size (nat \<circ> abs) (fst (snd p))  p2        ?(nat \<circ> abs) ...coeff?!?

   3) \<lambda>p. length (fst (snd p))                 p2

   4) \<lambda>p. size (snd p)                         (p2, m)

   5) \<lambda>p. list_size (nat \<circ> abs) (fst p)        p1        ?(nat \<circ> abs) ...coeff?!?

   6) \<lambda>p. length (fst p)                       p1

   7) size

 Result matrix:

     1  2  3  4  5  6  7 

 a:  <= ?  <  <= ?  ?  <=

 b:  <= <= <= <= ?  ?  <= *)

 declare mod_up_gcd.simps [simp del] --"HENSEL_lifting_up"

 value "mod_up_gcd [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] 7 = [2, 6, 0, 2, 6]" 

 value "mod_up_gcd [8, 28, 22, -11, -14, 1, 2] [2, 6, 0, 2, 6] 7 = [2, 6, 0, 2, 6]" 

 value "mod_up_gcd [20, 15, 8, 6] [8, -2, -2, 3] 2 = [0, 1]"

   value "[20, 15, 8, 6] %|% [0, 1] = False"

   value "[8, -2, -2, 3] %|% [0, 1] = False"

 (* analogous to Integer.gcds 

   gcds :: int list \<Rightarrow> int

   gcds ns = d

   assumes True

   yields THE d. ((\<forall>n. List.member ns n \<longrightarrow> d dvd n) \<and>

     (\<forall>d'. (\<forall>n. List.member ns n \<and> d' dvd n) \<longrightarrow> d'modp \<le> d)) *)

 fun gcds :: "int list \<Rightarrow> int" where "gcds ns = List.fold gcd ns (List.hd ns)"

 value "gcds [6, 9, 12] = 3"

 value "gcds [6, -9, 12] = 3"

 value "gcds [8, 12, 16] = 4"

 value "gcds [-8, 12, -16] = 4"

 (* prim_poly :: unipoly \<Rightarrow> unipoly

    prim_poly p = pp

    assumes True

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

 fun primitive_up :: "unipoly \<Rightarrow> unipoly" where

 "primitive_up [c] = (if c = 0 then [0] else [1])" |

 "primitive_up p =

   (let d = gcds p

   in

     if d = 1 then p else p div_ups d)"

 value "primitive_up [12, 16, 32, 44] = [3, 4, 8, 11]"

 value "primitive_up [4, 5, 12] =  [4, 5, 12]"

 value "primitive_up [0] = [0]"

 value "primitive_up [6] = [1]"

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

 (* 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: "p" is "prime", not "poly" ! *)

 function try_new_prime_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly" 

 where 

 "try_new_prime_up             a          b          d      M      P     g          p   = 

 (if P > M then g else 

   let p = p next_prime_not_dvd d;

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

                                 %* (int (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] else try_new_prime_up a b d M p g_p p

     else 

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

         let 

           P = P * p;

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

 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps"

 termination try_new_prime_up (*by lexicographic_order LOOPS +PROOF:? / by size_change LOOPS*)

 (*apply (relation "measures 

   [\<lambda>(a, b, d, M, P, g, p). ???,

    \<lambda>(a, b, d, M, P, g, p). ???,

    \<lambda>(a, b, d, M, P, g, p). ???]")

 Calls (manually):

   a) (a, b, d, M, P, g, p) ~> (a, b, d, M, p,         g_p {length g_p <= length a b}, p' {> p})

   b) (a, b, d, M, P, g, p) ~> (a, b, d, M, P,         g                             , p' {> p})

   c) (a, b, d, M, P, g, p) ~> (a, b, d, M, P {=P*p'}, g   {length g   <= length a b}, p' {> p})

 *)

 sorry

 (* 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: "p" is "prime", not "poly" ! *)

 function HENSEL_lifting_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly" where

 "HENSEL_lifting_up a b d M p = 

 (let

   p = p next_prime_not_dvd d;

   g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (int (d mod p))) mod_up p) p (*see above*)

 in 

   if deg_up g_p = 0 then [1] else 

     (let 

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

     in

       if (g %|% a) \<and> (g %|% b) then g else HENSEL_lifting_up a b d M p))"

 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps"

 termination HENSEL_lifting_up (*by lexicographic_order LOOPS +PROOF:?next_prime_not_dvd / by size_change LOOPS*)

 sorry 

 (*Calls (manually):

   a) (a, b, d, M, p) ~> (a, b, d, M, p {> p next_prime_not_dvd d})

 *)

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

 function gcd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" where

 "gcd_up a b = 

 (let d = \<bar>gcd (lcoeff_up a) (lcoeff_up b)\<bar>

 in

   if a = b then a else

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

 by pat_completeness auto --"simp del: lcoeff_up.simps ?+ others?"

 termination by lexicographic_order (*works*)

 value "gcd_up [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] = [-2, -1, 1]"

 (*     gcd    (-1 + x^2) (x + x^2) = (1 + x) ...*)

 value "gcd_up [-1, 0 ,1] [0, 1, 1] = [1, 1]"

 (*

 print_configs

 declare [[simp_trace_depth_limit = 99]]

 declare [[simp_trace = true]]

 using [[simp_trace_depth_limit = 99]]

 using [[simp_trace = true]]

 *)

 end