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

subsection {* a variant for div *}

fun 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 {* Meindl:  infix div' ... div' didnt work, infix didnt work *}

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

value " 5 div2  2"    -- "=  2"

value "-5 div2  2"    -- "= -2"

value "-5 div2 -2"    -- "=  2"

value " 5 div2 -2"    -- "= -2"

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

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

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

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

(* modi is just local for mod_inv *)

fun mod_inv :: "int => int => int"

  where "mod_inv r m = modi r r m 1"

text {*"----------- fun mod_inv --------------------------------"*}

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"

fun mod_div :: "int => int => int => int"

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

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

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

  by auto

  termination sorry

(* root1 is only local to approx_root *)

fun approx_root :: "int => int"

  where "approx_root a = root1 a 1"

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

text {*"----------- fun chinese_remainder ----------------------"*}

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 /\ n =< (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_prims is local to make_primes *)

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

  where "make_primes primes last_p n =

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

    then

      if is_prime primes (last_p + 2) 

      then make_primes (primes @ [(last_p + 2)]) (last_p + 2) n

      else make_primes  primes (last_p + 2) n

    else primes)"

  by pat_completeness (simp del: is_prime.simps)

  termination sorry

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

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

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

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

ML {* 

val ps = [2,3,5,7,11];

nth ps (List.length ps - 1);

*}

thm primes_upto.simps

thm primes_upto_def

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

text {*

  :

  by pat_completeness (simp del: is_prime.simps)

  by pat_completeness (simp del: primes_upto.simps)

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

  by pat_completeness (simp (no_asm))

  by simp

  by pat_completeness auto 

  by pat_completeness force 

  by pat_completeness blast

  :

function next_prime_not_dividing :: "int \<Rightarrow> int \<Rightarrow> int" (infixl "next'_prime'_not'_dividing" 70)

  where "p next_prime_not_dividing n = 

    (let

      ps = primes_upto (p + 1);

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

    in

      if n mod nxt \<noteq> 0

      then nxt

      else nxt next_prime_not_dividing n)"

  using [[simp_trace_depth_limit=999]]

  using [[simp_trace=true]]

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

  termination sorry

*}

end