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

 theory GCD_Poly_FP 

 imports "~~/src/HOL/Library/Polynomial" "~~/src/HOL/Number_Theory/Primes"

 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"

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

 subsection {* a variant for div *}

 (* 5 div 2 = 2; ~5 div 2 = ~3; but ~5 div2 2 = ~2; *)

 definition div2 :: "int => int => 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)"

 text {* div' didnt work *}

 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"

 subsection {* modulo calculations for integers *}

 (* modi is just local for mod_inv *)

 function modi :: "int => int => int => int => int" where

   "modi r rold m rinv = 

     (if rinv < m 

       then 

         if r mod m = 1 

         then rinv 

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

       else 0)"

   by auto

   termination sorry

 (* mod_inv :: int -> nat 

    mod_inv r m = s

    assumes 

    yields SOME s. r * s mod m = 1 *)

 fun mod_inv :: "int => int => int" 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"

 (* mod_div :: int -> int -> nat -> nat

    mod_div a b m = c

    assumes 

    yields SOME c. a * b^(-1) mod m = c *)

 definition mod_div :: "int => int => int => int" where

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

 (* root1 is just local to approx_root *)

 function root1 :: "int => int => int" where

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

 by auto

 termination sorry

 (* required just for one approximation:

    approx_root :: nat -> nat

    approx_root a = r

    assumes 

    yields SOME r. (r-1) * (r-1) < a  \<and>  a \<le> r * r*)

 definition approx_root :: "int => int" where

   "approx_root a = root1 a 1"

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

 definition chinese_remainder :: "int => int => int => int => int" where

   "chinese_remainder r1 m1 r2 m2 = 

     (r1 mod m1) + ((r2 - (r1 mod m1)) * (mod_inv m1 m2) mod m2) * m1"

 value "chinese_remainder 17 9 3 4  = 35"

 value "chinese_remainder  7 2 6 11 = 17"

 subsection {* operations with primes *}

 (* assumes: 'primes' contains all of first primes \<and> n \<le> (max primes)^2 *)

 fun is_prime :: "int list \<Rightarrow> int \<Rightarrow> bool" where

   "is_prime primes n = 

     (if List.length primes > 0

     then 

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

       then False

       else is_prime (List.tl primes) n

     else True)"

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

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

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

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

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

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

 (* make_primes is just local to primes_upto only *)

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

   "make_primes ps last_p n =

     (if (List.nth ps (List.length ps - 1)) < n

     then

       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

     else ps)"

 by pat_completeness (simp del: is_prime.simps)

 termination sorry

 (* primes_upto :: int \<Rightarrow> int list

    primes_upto n = [2, 3, 5, 7, .., p_k]                 -- sloppy formulations

      assumes 0 < n

      yields SOME k. ([2, 3, 5, 7, .., p_k] contains all primes \<le> p_k) \<and> p_(k-1) < n \<and> n \<le> p_k *)

 fun primes_upto :: "int \<Rightarrow> int list" where

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

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

 text {* FIXME next_prime_not_dvd loops with all the following:

 (* find a prime greater p not dividing the number n *)

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

  "p next_prime_not_dvd n = 

   (let

     ps = primes_upto (p + 1);

     nxt = List.nth ps (List.length ps - 1)

   in

     if n mod nxt \<noteq> 0

     then nxt

     else 77)"

 by pat_completeness (simp del: primes_upto.simps is_prime.simps) 

 termination sorry

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

   "p next_prime_not_dvd n = 

     (let

       ps = primes_upto (p + 1);

       nxt = List.nth ps (List.length ps - 1)

     in

       if n mod nxt \<noteq> 0

       then nxt

       else nxt next_prime_not_dvd n)"

 by pat_completeness (simp del: is_prime.simps make_primes.simps) 

 termination sorry

 ----

 loops                  : by pat_completeness auto

 loops                  : by pat_completeness simp

 Failed to finish proof : by pat_completeness (simp del: is_prime.simps make_primes.simps)

 Failed to finish proof : by pat_completeness (simp del: make_primes.simps)

 *}

 declare [[simp_trace_depth_limit=99]]

 declare [[simp_trace=true]]

 ML {*11111*}

 thm make_primes.simps

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

   "p next_prime_not_dvd n = 

     (let

       ps = primes_upto (p + 1);

       nxt = List.nth ps (List.length ps - 1)

     in

       if n mod nxt \<noteq> 0

       then nxt

       else nxt next_prime_not_dvd n)"

 by pat_completeness (simp del: make_primes.simps)

 (*---------------------------------------------------------------------------------------------

 subsection {* auxiliary functions 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 *)

 function lc :: "unipoly \<Rightarrow> int" where

   "lc uvp =

     (if List.nth uvp (List.length uvp - 1) \<noteq> 0

     then List.nth uvp (List.length uvp - 1)

     else lc (nth_drop (List.length uvp - 1) uvp))"

 by auto 

 termination sorry

 text {*"----------- fun lc -------------------------------------"*}

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

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

 (* degree *)

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

   "deg uvp =

     (if nth uvp (List.length uvp - 1) \<noteq> 0

     then List.length uvp - 1

     else deg (nth_drop (List.length uvp - 1) uvp))"

 by auto 

 termination sorry

 (* TODO: ?let l = length uvp - 1? more efficient??? compare drop_lc0 !!! *)

 text {*"----------- fun deg ------------------------------------"*}

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

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

 (* drop leading coefficients equal 0 *)

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

   "drop_lc0 [] = []"

 | "drop_lc0 uvp = 

     (let l = List.length uvp - 1

     in

       if nth uvp l = 0

       then drop_lc0 (nth_drop l uvp)

       else uvp)"

 text {*"----------- fun drop_lc0 -------------------------------";*}

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

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

 (* norm is just local to normalise *)

 fun norm :: "unipoly \<Rightarrow> unipoly \<Rightarrow> int (*nat?*) \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> unipoly" where

   "norm pol nrm m lcp i = 

     (if i = 0

         then [mod_div (nth pol i) lcp m] @ nrm

         else norm pol ([mod_div (nth pol i) lcp m] @  nrm) m lcp (i - 1))"

 (* normalise a unipoly such that the lc 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> int (*nat?*) \<Rightarrow> unipoly" where

   "normalise p m = 

     (let

      p = drop_lc0 p;

      lcp = lc p

     in 

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

 text {*"----------- fun normalise ------------------------------";*}

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

 (* scalar multiplication, %* in GCD_Poly.thy *)

 fun  mult_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "mult'_ups" 70) where

  "p mult_ups m = map (op * m) p"

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

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

 (* scalar divison, %/ in GCD_Poly.thy *)

 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" 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,  %-% in GCD_Poly.thy *)

 definition minus_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" (infixl "minus'_up" 70) where

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

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

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

 (* dvd for univariate polynomials *)

 function dvd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> bool" (infixl "dvd'_up" 70) where

 (* FIXME this makes induction proof loop .......

   "[d] dvd_up [i] =

     (if (\<bar>d\<bar> \<le> \<bar>i\<bar>) & (i mod d = 0) then True else False)"

 |*) "d dvd_up p =

     (let 

       d = drop_lc0 d;

       p = drop_lc0 p;

       d000 = (List.replicate (List.length p - List.length d) 0) @ d;

       quot = (lc p) div2 (lc d000);

       rest = drop_lc0 (p minus_up (d000 div_ups quot))

     in 

       if rest = [] then True

       else 

         (if quot \<noteq> 0 & List.length d \<le> List.length rest

         then d dvd_up rest

         else False))"

 by pat_completeness auto 

 termination sorry

 thm dvd_up.simps

 thm dvd_up.induct

 (* 

 value "[2, 3] dvd_up [8, 22, 15] = True"

 ... FIXME loops, incorrect: ALL SHOULD EVALUATE TO TRUE .......

 *)

 value "[4] dvd_up [6]                 = False"

 value "[8] dvd_up [16, 0]             = True"

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

 value "[8, 0] dvd_up [16]             = True"

 (* centr is just local to poly_centr *)

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

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

 (* TODO *)

 definition poly_centr :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" where

   "poly_centr p m =

    (let

       mi = m div2 2;

       mid = if m mod mi = 0 then mi else mi + 1

    in map (centr m mid) p)"

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

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

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

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

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

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

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

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

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

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

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

 (* TODO *)

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

   "sum_lmb p exp = 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"

 (* TODO *)

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

   "landau_mignotte_bound p1 p2 = 

     ((power 2 (min (deg p1) (deg p2))) * (nat (gcd (lc p1) (lc p2))) * 

     (nat (min (abs ((approx_root (sum_lmb p1 2)) div2 -(lc p1)))

                 (abs ((approx_root (sum_lmb p2 2)) div2 -(lc 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_poly *)

 fun pair :: "unipoly \<Rightarrow> unipoly \<Rightarrow> ((int * 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)"

 (* chinese_rem is just local to chinese_remainder_poly *)

 fun chinese_rem :: "int * int \<Rightarrow> int * int \<Rightarrow> int" where

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

 (* TODO *)

 fun chinese_remainder_poly :: "int * int \<Rightarrow> unipoly * unipoly \<Rightarrow> unipoly" where

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

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

 (* mod_pol is just local to mod_poly *)

 function mod_pol :: "unipoly \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> unipoly" where

   "mod_pol p m n = (if n < List.length p

       then mod_pol ((List.tl p) @ [(List.hd p) mod m]) m (n + 1)

       else p)"

 by pat_completeness auto 

 termination sorry

 (* TODO *)

 definition mod_poly :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "mod'_poly" 70) where

   "p mod_poly m = mod_pol p m 0"

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

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

 (* euclidean algoritm in Z_p[x].

    mod_poly_gcd :: unipoly -> unipoly -> int -> unipoly

    mod_poly_gcd p1 p2 m = g

      assumes 

      yields \<exists> g. gcd((p1 mod m),(p2 mod m)) = g *)

 function mod_poly_gcd :: "unipoly \<Rightarrow> unipoly  \<Rightarrow> int \<Rightarrow> unipoly" where

  "mod_poly_gcd p1 p2 m = 

    (let 

       p1m = p1 mod_poly m;

       p2m = drop_lc0 (p2 mod_poly m);

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

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

       rest = drop_lc0 ((p1m minus_up (p2n mult_ups quot)) mod_poly m)

     in 

    if rest = [] 

         then p2

         else

           if length rest < List.length p2

           then mod_poly_gcd p2 rest m 

           else mod_poly_gcd rest p2 m)"

 by pat_completeness auto 

 termination sorry

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

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

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

 (*----------------------------------------------------------------------------*)

 value "xxx" -- "= "

 print_configs

 declare [[simp_trace_depth_limit=99]]

 declare [[simp_trace=true]]

 ---------------------------------------------------------------------------------------------*)

 end