Theory Prime_Factorization

(*  
    Author:      René Thiemann 
                 Akihisa Yamada
    License:     BSD
*)
section ‹Prime Factorization›

text ‹This theory contains not-completely naive algorithms to test primality
  and to perform prime factorization. More precisely, it corresponds to prime factorization
  algorithm A in Knuth's textbook cite"Knuth".›

theory Prime_Factorization
imports
  "HOL-Computational_Algebra.Primes"
  Missing_List
  Missing_Multiset
begin

subsection ‹Definitions›

definition primes_1000 :: "nat list" where
  "primes_1000 = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101,
  103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
  211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
  331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
  449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577,
  587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
  709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839,
  853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983,
  991, 997]"

lemma primes_1000: "primes_1000 = filter prime [0..<1001]"
  by eval

definition next_candidates :: "nat  nat × nat list" where
  "next_candidates n = (if n = 0 then (1001,primes_1000) else (n + 30, 
    [n,n+2,n+6,n+8,n+12,n+18,n+20,n+26]))"

definition "candidate_invariant n = (n = 0  n mod 30 = (11 :: nat))"

partial_function (tailrec) remove_prime_factor :: "nat  nat  nat list  nat × nat list " where
  [code]: "remove_prime_factor p n ps = (case Euclidean_Rings.divmod_nat n p of (n',m)  
     if m = 0 then remove_prime_factor p n' (p # ps) else (n,ps))" 
  
partial_function (tailrec) prime_factorization_nat_main 
  :: "nat  nat  nat list  nat list  nat list" where
  [code]: "prime_factorization_nat_main n j is ps = (case is of 
     []  
       (case next_candidates j of (j,is)  prime_factorization_nat_main n j is ps)
   | (i # is)  (case Euclidean_Rings.divmod_nat n i of (n',m)  
       if m = 0 then case remove_prime_factor i n' (i # ps)
       of (n',ps')  if n' = 1 then ps' else 
         prime_factorization_nat_main n' j is ps'
       else if i * i  n then prime_factorization_nat_main n j is ps
       else (n # ps)))"

partial_function (tailrec) prime_nat_main 
  :: "nat  nat  nat list  bool" where
  [code]: "prime_nat_main n j is = (case is of 
    []  (case next_candidates j of (j,is)  prime_nat_main n j is)
  | (i # is)  (if i dvd n then i  n else if i * i  n then prime_nat_main n j is
    else True))"

definition prime_nat :: "nat  bool" where
  "prime_nat n  if n < 2 then False else ― ‹TODO: integrate precomputed map›
     case next_candidates 0 of (j,is)  prime_nat_main n j is"

definition prime_factorization_nat :: "nat  nat list" where
  "prime_factorization_nat n  rev (if n < 2 then [] else 
    case next_candidates 0 of (j,is)  prime_factorization_nat_main n j is [])"

definition divisors_nat :: "nat  nat list" where 
  "divisors_nat n  if n = 0 then [] else 
     remdups_adj (sort (map prod_list (subseqs (prime_factorization_nat n))))"

definition divisors_int_pos :: "int  int list" where
  "divisors_int_pos x  map int (divisors_nat (nat (abs x)))"

definition divisors_int :: "int  int list" where
  "divisors_int x  let xs = divisors_int_pos x in xs @ (map uminus xs)"

subsection ‹Proofs›

lemma remove_prime_factor: assumes res: "remove_prime_factor i n ps = (m,qs)"
  and i: "i > 1"
  and n: "n  0"
  shows " rs. qs = rs @ ps  n = m * prod_list rs  ¬ i dvd m  set rs  {i}"
  using res n
proof (induct n arbitrary: ps rule: less_induct)
  case (less n ps)
  obtain n' mo where dm: "Euclidean_Rings.divmod_nat n i = (n',mo)" by force
  hence n': "n' = n div i" and mo: "mo = n mod i" by (auto simp: Euclidean_Rings.divmod_nat_def)
  from less(2)[unfolded remove_prime_factor.simps[of i n] dm]
  have res: "(if mo = 0 then remove_prime_factor i n' (i # ps) else (n, ps)) = (m, qs)" by auto
  from less(3) have n: "n  0" by auto
  with n' i have "n' < n" by auto
  note IH = less(1)[OF this]
  show ?case
  proof (cases "mo = 0")
    case True
    with mo n' have n: "n = n' * i" by auto
    with n  0 have n': "n'  0" by auto
    from True res have "remove_prime_factor i n' (i # ps) = (m,qs)" by auto
    from IH[OF this n'] obtain rs where 
      "qs = rs @ i # ps" and "n' = m * prod_list rs  ¬ i dvd m  set rs  {i}" by auto
    thus ?thesis
      by (intro exI[of _ "rs @ [i]"], unfold n, auto)
  next
    case False
    with mo have i_n: "¬ i dvd n" by auto
    from False res have id: "m = n" "qs = ps" by auto
    show ?thesis unfolding id using i_n by auto
  qed
qed 

lemma prime_sqrtI: assumes n: "n  2" 
  and small: " j. 2  j  j < i  ¬ j dvd n"
  and i: "¬ i * i  n"
  shows "prime (n::nat)" unfolding prime_nat_iff
proof (intro conjI impI allI)
  show "1 < n" using n by auto
  fix j
  assume jn: "j dvd n" 
  from jn obtain k where njk: "n = j * k" unfolding dvd_def by auto
  with 1 < n have jn: "j  n" by (metis dvd_imp_le jn neq0_conv not_less0)
  show "j = 1  j = n"
  proof (rule ccontr)
    assume "¬ ?thesis"
    with njk n have j1: "j > 1  j  n" by simp
    have " j k. 1 < j  j  k  n = j * k"
    proof (cases "j  k")
      case True
      thus ?thesis unfolding njk using j1 by blast
    next
      case False
      show ?thesis by (rule exI[of _ k], rule exI[of _ j], insert 1 < n j1 njk False, auto)
        (metis Suc_lessI mult_0_right neq0_conv)
    qed
    then obtain j k where j1: "1 < j" and jk: "j  k" and njk: "n = j * k" by auto
    with small[of j] have ji: "j  i" unfolding dvd_def by force
    from mult_mono[OF ji ji] have "i * i  j * j" by auto
    with i have "j * j > n" by auto
    from this[unfolded njk] have "k < j" by auto
    with jk show False by auto
  qed
qed

lemma candidate_invariant_0: "candidate_invariant 0"
  unfolding candidate_invariant_def by auto

lemma next_candidates: assumes res: "next_candidates n = (m,ps)"
  and n: "candidate_invariant n"
  shows "candidate_invariant m" "sorted ps" "{i. prime i  n  i  i < m}  set ps" 
    "set ps  {2..}  {n..<m}" "distinct ps" "ps  []" "n < m"
  unfolding candidate_invariant_def
proof -
  note res = res[unfolded next_candidates_def]
  note n = n[unfolded candidate_invariant_def]
  show "m = 0  m mod 30 = 11" using res n by (auto split: if_splits)
  show "sorted ps" using res n by (auto split: if_splits simp: primes_1000_def sorted2_simps simp del: sorted_wrt.simps(2))
  show "set ps  {2..}  {n..<m}" using res n by (auto split: if_splits simp: primes_1000_def)
  show "distinct ps" using res n by (auto split: if_splits simp: primes_1000_def)
  show "ps  []" using res n by (auto split: if_splits simp: primes_1000_def)
  show "n < m" using res by (auto split: if_splits)
  show "{i. prime i  n  i  i < m}  set ps"
  proof (cases "n = 0")
    case True
    hence *: "m = 1001" "ps = primes_1000" using res by auto
    show ?thesis unfolding * True primes_1000 by auto
  next
    case False
    hence n: "n mod 30 = 11" and m: "m = n + 30" and ps: "ps = [n,n+2,n+6,n+8,n+12,n+18,n+20,n+26]" 
      using res n by auto
    {
      fix i
      assume *: "prime i" "n  i" "i < n + 30" "i  set ps"
      from n * have i11: "i  11" by auto
      define j where "j = i - n" 
      have i: "i = n + j" using n  i j_def by auto
      have "i mod 30 = (j + n) mod 30" using n  i unfolding j_def by simp
      also have " = (j mod 30 + n mod 30) mod 30"
        by (simp add: mod_simps)
      also have " = (j mod 30 + 11) mod 30" unfolding n by simp
      finally have i30: "i mod 30 = (j mod 30 + 11) mod 30" by simp
      have 2: "2 dvd (30 :: nat)" and 112: "11 mod (2 :: nat) = 1" by simp_all
      have "(j + 11) mod 2 = (j + 1) mod 2"
        by (rule mod_add_cong) simp_all
      with arg_cong [OF i30, of "λj. j mod 2"]
      have 2: "i mod 2 = (j mod 2 + 1) mod 2"
        by (simp add: mod_simps mod_mod_cancel [OF 2])
      have 3: "3 dvd (30 :: nat)" and 113: "11 mod (3 :: nat) = 2" by simp_all
      have "(j + 11) mod 3 = (j + 2) mod 3"
        by (rule mod_add_cong) simp_all
      with arg_cong [OF i30, of "λ j. j mod 3"] have 3: "i mod 3 = (j mod 3 + 2) mod 3"
        by (simp add: mod_simps mod_mod_cancel [OF 3])
      have 5: "5 dvd (30 :: nat)" and 115: "11 mod (5 :: nat) = 1" by simp_all
      have "(j + 11) mod 5 = (j + 1) mod 5"
        by (rule mod_add_cong) simp_all
      with arg_cong [OF i30, of "λ j. j mod 5"] have 5: "i mod 5 = (j mod 5 + 1) mod 5"
        by (simp add: mod_simps mod_mod_cancel [OF 5])
      
      from n *(2-)[unfolded ps i, simplified] have 
        "j  {1,3,5,7,9,11,13,15,17,19,21,23,25,27,29}  j  {4,10,16,22,28}  j  {14,24}"
        (is "j  ?j2  j  ?j3  j  ?j5")
        by simp presburger
      moreover
      {
        assume "j  ?j2"
        hence "j mod 2 = 1" by auto
        with 2 have "i mod 2 = 0" by auto
        with i11 have "2 dvd i" "i  2" by auto 
        with *(1) have False unfolding prime_nat_iff by auto
      }
      moreover
      {
        assume "j  ?j3"
        hence "j mod 3 = 1" by auto
        with 3 have "i mod 3 = 0" by auto
        with i11 have "3 dvd i" "i  3" by auto 
        with *(1) have False unfolding prime_nat_iff by auto
      }
      moreover
      {
        assume "j  ?j5"
        hence "j mod 5 = 4" by auto
        with 5 have "i mod 5 = 0" by auto
        with i11 have "5 dvd i" "i  5" by auto 
        with *(1) have False unfolding prime_nat_iff by auto
      }
      ultimately have False by blast
    }
    thus ?thesis unfolding m ps by auto
  qed
qed


lemma prime_test_iterate2: assumes small: " j. 2  j  j < (i :: nat)  ¬ j dvd n"
  and odd: "odd n"
  and n: "n  3"
  and i: "i  3" "odd i"
  and mod: "¬ i dvd n"
  and j: "2  j" "j < i + 2"
  shows "¬ j dvd n"
proof 
  assume dvd: "j dvd n"
  with small[OF j(1)] have "j  i" by linarith
  with dvd mod have "j > i" by (cases "i = j", auto)
  with j have "j = Suc i" by simp
  with i have "even j" by auto
  with dvd j(1) have "2 dvd n" by (metis dvd_trans)
  with odd show False by auto
qed

lemma prime_divisor: assumes "j  2" and "j dvd n" shows
  " p :: nat. prime p  p dvd j  p dvd n"
proof -
  let ?pf = "prime_factors j"
  from assms have "j > 0" by auto
  from prime_factorization_nat[OF this]
  have "j = (p?pf. p ^ multiplicity p j)" by auto
  with j  2 have "?pf  {}" by auto
  then obtain p where p: "p  ?pf" by auto  
  hence pr: "prime p" by auto
  define rem where "rem = (p?pf - {p}. p ^ multiplicity p j)"
  from p have mult: "multiplicity p j  0"
    by (auto simp: prime_factors_multiplicity)
  have "finite ?pf" by simp
  have "j = (p?pf. p ^ multiplicity p j)" by fact
  also have "?pf = (insert p (?pf - {p}))" using p by auto
  also have "(pinsert p (?pf - {p}). p ^ multiplicity p j) = 
    p ^ multiplicity p j * rem" unfolding rem_def
    by (subst prod.insert, auto)
  also have " = p * (p ^ (multiplicity p j - 1) * rem)" using mult 
    by (cases "multiplicity p j", auto)
  finally have pj: "p dvd j" unfolding dvd_def by blast
  with j dvd n have "p dvd n" by (metis dvd_trans)
  with pj pr show ?thesis by blast
qed

lemma prime_nat_main: "ni = (n,i,is)  i  2  n  2 
  ( j. 2  j  j < i  ¬ (j dvd n)) 
  ( j. i  j  j < jj  prime j  j  set is)  i  jj 
  sorted is  distinct is  candidate_invariant jj  set is  {i..<jj}  
  res = prime_nat_main n jj is  
  res = prime n"
proof (induct ni arbitrary: n i "is" jj res rule: wf_induct[OF 
    wf_measures[of "[λ (n,i,is). n - i, λ (n,i,is). if is = [] then 1 else 0]"]])
  case (1 ni n i "is" jj res)
  note res = 1(12)
  from 1(3-4) have i: "i  2" and i2: "Suc i  2" and n: "n  2" by auto
  from 1(5) have dvd: " j. 2  j  j < i  ¬ j dvd n" .
  from 1(7) have ijj: "i  jj" .
  note sort_dist = 1(8-9)
  have "is": " j. i  j  j < jj  prime j  j  set is" by (rule 1(6))
  note simps = prime_nat_main.simps[of n jj "is"]
  note IH = 1(1)[rule_format, unfolded 1(2), OF _ refl]
  show ?case
  proof (cases "is")
    case Nil
    obtain jjj iis where can: "next_candidates jj = (jjj,iis)" by force
    from res[unfolded simps, unfolded Nil can split] have res: "res = prime_nat_main n jjj iis" by auto
    from next_candidates[OF can 1(10)] have can: 
      "sorted iis" "distinct iis" "candidate_invariant jjj" 
      "{i. prime i  jj  i  i < jjj}  set iis" "set iis  {2..}  {jj..<jjj}"
      "iis  []" "jj < jjj" by blast+
    from can ijj have "i  jjj" by auto
    note IH = IH[OF _ i n dvd _ this can(1-3) _ res]
    show ?thesis
    proof (rule IH, force simp: Nil can(6))
      fix x
      assume ix: "i  x" and xj: "x < jjj" and px: "prime x"
      from "is"[OF ix _ px] Nil have jx: "jj  x" by force
      with can(4) xj px show "x  set iis" by auto
    qed (insert can(5) ijj, auto)
  next
    case (Cons i' iis)
    with res[unfolded simps]
    have res: "res = (if i' dvd n then n  i' else if i' * i'  n then prime_nat_main n jj iis else True)" 
      by simp
    from 1(11) Cons have iis: "set iis  {i..<jj}" and i': "i  i'" "i' < jj" "Suc i'  jj" by auto
    from sort_dist have sd_iis: "sorted iis" "distinct iis" and "i'  set iis" by(auto simp: Cons)
    from sort_dist(1) have "set iis  {i'..}" by(auto simp: Cons)
    with iis have "set iis  {i'..<jj}" by force
    with i'  set iis  have iis: "set iis  {Suc i'..<jj}"
      by (auto, case_tac "x = i'", auto)
    {
      fix j
      assume j: "2  j" "j < i'"
      have "¬ j dvd n"
      proof
        assume "j dvd n"
        from prime_divisor[OF j(1) this] obtain p where 
          p: "prime p" "p dvd j" "p dvd n" by auto
        have pj: "p  j"
          by (rule dvd_imp_le[OF p(2)], insert j, auto)
        have p2: "2  p" using p(1) by (rule prime_ge_2_nat)
        from dvd[OF p2] p(3) have pi: "p  i" by force
        from pj j(2) i' "is"[OF pi _ p(1)] have "p  set is" by auto
        with sorted is have "i'  p" by(auto simp: Cons)
        with pj j(2) show False by arith
      qed
    } note dvd = this
    from i' i have i'2: "2  Suc i'" by auto
    note IH = IH[OF _ i'2 n _ _ i'(3) sd_iis 1(10) iis]
    show ?thesis
    proof (cases "i' dvd n")
      case False note dvdi = this
      {
        fix j
        assume j: "2  j" "j < Suc i'"
        have "¬ j dvd n"
        proof (cases "j = i'")
          case False
          with j have "j < i'" by auto
          from dvd[OF j(1) this] show ?thesis .
        qed (insert False, auto)
      } note dvds = this
      show ?thesis
      proof (cases "i' * i'  n")
        case True note iin = this
        with res False have res: "res = prime_nat_main n jj iis" by auto
        from iin have i_n: "i' < n"
          using dvd dvdi n nat_neq_iff dvd_refl by blast
        {
          fix x
          assume "Suc i'  x" "x < jj" "prime x"
          hence "i  x" "x < jj" "prime x" using i' by auto
          from "is"[OF this] have "x  set is" . 
          with Suc i'  x have "x  set iis" unfolding Cons by auto
        } note iis = this
        show ?thesis 
          by (rule IH[OF _ dvds iis res], insert i_n i', auto)
      next
        case False
        with res dvdi have res: "res = True" by auto
        have n: "prime n" 
          by (rule prime_sqrtI[OF n dvd False])
        thus ?thesis unfolding res by auto
      qed
    next
      case True
      have "i'  2" using i i' by auto
      from i' dvd n obtain k where "n = i' * k" ..
      with n have "k  0" by (cases "k = 0", auto)
      with n = i' * k have *: "i' < n  i' = n"
        by auto
      with True res have "res  i' = n"
        by auto
      also have " = prime n"
      using * proof
        assume "i' < n"
        with i'  2 i' dvd n have "¬ prime n"
          by (auto simp add: prime_nat_iff)
        with i' < n show ?thesis
          by auto
      next
        assume "i' = n"
        with dvd n have "prime n"
          by (simp add: prime_nat_iff')
        with i' = n show ?thesis 
          by auto
      qed
      finally show ?thesis .
    qed
  qed
qed

lemma prime_factorization_nat_main: "ni = (n,i,is)  i  2  n  2 
  ( j. 2  j  j < i  ¬ (j dvd n))  
  ( j. i  j  j < jj  prime j  j  set is)  i  jj 
  sorted is  distinct is  candidate_invariant jj  set is  {i..<jj}  
  res = prime_factorization_nat_main n jj is ps  
   qs. res = qs @ ps  Ball (set qs) prime  n = prod_list qs"
proof (induct ni arbitrary: n i "is" jj res ps rule: wf_induct[OF 
  wf_measures[of "[λ (n,i,is). n - i, λ (n,i,is). if is = [] then 1 else 0]"]])
  case (1 ni n i "is" jj res ps)
  note res = 1(12)
  from 1(3-4) have i: "i  2" and i2: "Suc i  2" and n: "n  2" by auto
  from 1(5) have dvd: " j. 2  j  j < i  ¬ j dvd n" .
  from 1(7) have ijj: "i  jj" .
  note sort_dist = 1(8-9)
  have "is": " j. i  j  j < jj  prime j  j  set is" by (rule 1(6))
  note simps = prime_factorization_nat_main.simps[of n jj "is"]
  note IH = 1(1)[rule_format, unfolded 1(2), OF _ refl]
  show ?case
  proof (cases "is")
    case Nil
    obtain jjj iis where can: "next_candidates jj = (jjj,iis)" by force
    from res[unfolded simps, unfolded Nil can split] have res: "res = prime_factorization_nat_main n jjj iis ps" by auto
    from next_candidates[OF can 1(10)] have can: 
      "sorted iis" "distinct iis" "candidate_invariant jjj" 
      "{i. prime i  jj  i  i < jjj}  set iis" "set iis  {2..}  {jj..<jjj}"
      "iis  []" "jj < jjj" by blast+
    from can ijj have "i  jjj" by auto
    note IH = IH[OF _ i n dvd _ this can(1-3) _ res]
    show ?thesis
    proof (rule IH, force simp: Nil can(6))
      fix x
      assume ix: "i  x" and xj: "x < jjj" and px: "prime x"
      from "is"[OF ix _ px] Nil have jx: "jj  x" by force
      with can(4) xj px show "x  set iis" by auto
    qed (insert can(5) ijj, auto)
  next
    case (Cons i' iis)
    obtain n' m where dm: "Euclidean_Rings.divmod_nat n i' = (n',m)" by force
    hence n': "n' = n div i'" and m: "m = n mod i'" by (auto simp: Euclidean_Rings.divmod_nat_def)
    have m: "(m = 0) = (i' dvd n)" unfolding m by auto
    from Cons res[unfolded simps] dm m n'
    have res: "res = (
       if i' dvd n then case remove_prime_factor i' (n div i') (i' # ps) of
            (n', ps')  if n' = 1 then ps' else prime_factorization_nat_main n' jj iis ps'
       else if i' * i'  n then prime_factorization_nat_main n jj iis ps else n # ps)" 
      by simp
    from 1(11) i Cons have iis: "set iis  {i..<jj}" and i': "i  i'" "i' < jj" "Suc i'  jj" "i' > 1" by auto
    from sort_dist have sd_iis: "sorted iis" "distinct iis" and "i'  set iis" by(auto simp: Cons)
    from sort_dist(1) Cons have "set iis  {i'..}" by(auto)
    with iis have "set iis  {i'..<jj}" by force
    with i'  set iis  have iis: "set iis  {Suc i'..<jj}"
      by (auto, case_tac "x = i'", auto)
    {
      fix j
      assume j: "2  j" "j < i'"
      have "¬ j dvd n"
      proof
        assume "j dvd n"
        from prime_divisor[OF j(1) this] obtain p where 
          p: "prime p" "p dvd j" "p dvd n" by auto
        have pj: "p  j"
          by (rule dvd_imp_le[OF p(2)], insert j, auto)
        have p2: "2  p" using p(1) by (rule prime_ge_2_nat)
        from dvd[OF p2] p(3) have pi: "p  i" by force
        from pj j(2) i' "is"[OF pi _ p(1)] have "p  set is" by auto
        with sorted is have "i'  p" by (auto simp: Cons)
        with pj j(2) show False by arith
      qed
    } note dvd = this
    from i' i have i'2: "2  Suc i'" by auto
    note IH = IH[OF _ i'2 _ _ _ i'(3) sd_iis 1(10) iis]
    {
      fix x
      assume "Suc i'  x" "x < jj" "prime x"
      hence "i  x" "x < jj" "prime x" using i' by auto
      from "is"[OF this] have "x  set is" . 
      with Suc i'  x have "x  set iis" unfolding Cons by auto
    } note iis = this
    show ?thesis
    proof (cases "i' dvd n")
      case False note dvdi = this
      {
        fix j
        assume j: "2  j" "j < Suc i'"
        have "¬ j dvd n"
        proof (cases "j = i'")
          case False
          with j have "j < i'" by auto
          from dvd[OF j(1) this] show ?thesis .
        qed (insert False, auto)
      } note dvds = this
      show ?thesis
      proof (cases "i' * i'  n")
        case True note iin = this
        with res False have res: "res = prime_factorization_nat_main n jj iis ps" by auto
        from iin have i_n: "i' < n" using dvd dvdi n nat_neq_iff dvd_refl by blast 
        show ?thesis 
          by (rule IH[OF _ n dvds iis res], insert i_n i', auto)
      next
        case False
        with res dvdi have res: "res = n # ps" by auto
        have n: "prime n" 
          by (rule prime_sqrtI[OF n dvd False])
        thus ?thesis unfolding res by auto
      qed    
    next
      case True note i_n = this
      obtain n'' qs where rp: "remove_prime_factor i' (n div i') (i' # ps) = (n'',qs)" by force
      with res True 
      have res: "res = (if n'' = 1 then qs else prime_factorization_nat_main n'' jj iis qs)" by auto
      have pi: "prime i'" unfolding prime_nat_iff
      proof (intro conjI allI impI)
        show "1 < i'" using i' i by auto
        fix j
        assume ji: "j dvd i'"
        with i' i have j0: "j  0" by (cases "j = 0", auto)
        from ji i_n have jn: "j dvd n" by (metis dvd_trans)
        with dvd[of j] have j: "2 > j  j  i'" by linarith
        from ji 1 < i' have "j  i'" unfolding dvd_def 
          by (simp add: dvd_imp_le ji)
        with j j0 show "j = 1  j = i'" by linarith
      qed
      from True n' have id: "n = n' * i'" by auto
      from n id have "n'  0" by (cases "n = 0", auto)
      with id have "i'  n" by auto
      from remove_prime_factor[OF rp[folded n'] 1 < i' n'  0] obtain rs
        where qs: "qs = rs @ i' # ps" and n': "n' = n'' * prod_list rs" and i_n'': "¬ i' dvd n''" 
        and rs: "set rs  {i'}" by auto
      {
        fix j
        assume "j dvd n''"
        hence "j dvd n" unfolding id n' by auto
      } note dvd' = this
      show ?thesis
      proof (cases "n'' = 1")
        case False
        with res have res: "res = prime_factorization_nat_main n'' jj iis qs" 
          by simp
        from i i' have "i'  2" by simp
        from False n' n'  0 have n2: "n''  2" by (cases "n'' = 0"; auto)
        have lrs: "prod_list rs  0" using n' n'  0 by (cases "prod_list rs = 0", auto)
        with i'  2 have "prod_list rs * i'  2" by (cases "prod_list rs", auto)
        hence nn'': "n > n''" unfolding id n' using n2 by simp
        have "i'  n" unfolding id n' using pi False by fastforce
        with i'  n i' have "n > i" by auto
        with nn'' i i' have less: "n - i > n'' - Suc i'" by simp
        {
          fix j
          assume 2: "2  j" and ji: "j < Suc i'" 
          have "¬ j dvd n''" 
          proof (cases "j = i'")
            case False
            with ji have "j < i'" by auto
            from dvd' dvd[OF 2 this] show ?thesis by blast
          qed (insert i_n'', auto)
        }
        from IH[OF _ n2 this iis res] less obtain ss where 
          res: "res = ss @ qs  Ball (set ss) prime  n'' = prod_list ss" by auto
        thus ?thesis unfolding id n' qs using pi rs by auto
      next
        case True
        with res have res: "res = qs" by auto
        show ?thesis unfolding id n' res qs True using rs prime i'
          by (intro exI[of _ "rs @ [i']"], auto)
      qed
    qed
  qed
qed

lemma prime_nat[simp]: "prime_nat n = prime n"
proof (cases "n < 2")
  case True
  thus ?thesis unfolding prime_nat_def prime_nat_iff by auto
next
  case False
  hence n: "n  2" by auto
  obtain jj "is" where can: "next_candidates 0 = (jj,is)" by force
  from next_candidates[OF this candidate_invariant_0]
  have cann: "sorted is" "distinct is" "candidate_invariant jj" 
    "{i. prime i  0  i  i < jj}  set is"
    "set is  {2..}  {0..<jj}" "distinct is" "is  []" by auto  
  from cann have sub: "set is  {2..<jj}" by force
  with is  [] have jj: "jj  2" by (cases "is", auto)
  from n can have res: "prime_nat n = prime_nat_main n jj is"
    unfolding prime_nat_def by auto
  show ?thesis using prime_nat_main[OF refl le_refl n _ _ jj cann(1-3) sub res] cann(4) by auto
qed

lemma prime_factorization_nat: fixes n :: nat
  defines "pf  prime_factorization_nat n"
  shows "Ball (set pf) prime"
  and "n  0  prod_list pf = n"
  and "n = 0  pf = []"
proof -
  note pf = pf_def[unfolded prime_factorization_nat_def]
  have "Ball (set pf) prime  (n  0  prod_list pf = n)  (n = 0  pf = [])"
  proof (cases "n < 2")
    case True
    thus ?thesis using pf by auto
  next
    case False
    hence n: "n  2" by auto
    obtain jj "is" where can: "next_candidates 0 = (jj,is)" by force
    from next_candidates[OF this candidate_invariant_0]
    have cann: "sorted is" "distinct is" "candidate_invariant jj" 
      "{i. prime i  0  i  i < jj}  set is"
      "set is  {2..}  {0..<jj}" "distinct is" "is  []" by auto  
    from cann have sub: "set is  {2..<jj}" by force
    with is  [] have jj: "jj  2" by (cases "is", auto)
    let ?pfm = "prime_factorization_nat_main n jj is []"
    from pf[unfolded can] False
    have res: "pf = rev ?pfm" by simp
    from prime_factorization_nat_main[OF refl le_refl n _ _ jj cann(1-3) sub refl, of Nil] cann(4)
    have "Ball (set ?pfm) prime" "n = prod_list ?pfm" by auto
    thus ?thesis unfolding res using n by auto
  qed
  thus "Ball (set pf) prime" "n  0  prod_list pf = n" "n = 0  pf = []" by auto
qed

lemma prod_mset_multiset_prime_factorization_nat [simp]: 
  "(x::nat)  0  prod_mset (prime_factorization x) = x"
  by simp

(* TODO Move *)
lemma prime_factorization_unique'':
  fixes A :: "'a :: {factorial_semiring_multiplicative} multiset"
  assumes "p. p ∈# A  prime p"
  assumes "prod_mset A = normalize x"
  shows   "prime_factorization x = A"
proof -
  have "prod_mset A  0" by (auto dest: assms(1))
  with assms(2) have "x  0" by simp
  hence "prod_mset (prime_factorization x) = prod_mset A" 
    by (simp add: assms prod_mset_prime_factorization)
  with assms show ?thesis 
    by (intro prime_factorization_unique') auto
qed

lemma multiset_prime_factorization_nat_correct:
  "prime_factorization n = mset (prime_factorization_nat n)"
proof -
  note pf = prime_factorization_nat[of n]
  show ?thesis
  proof (cases "n = 0")
    case True
    thus ?thesis using pf(3) by simp
  next
    case False
    note pf = pf(1) pf(2)[OF False]
    show ?thesis
    proof (rule prime_factorization_unique'')
      show "prime p" if "p ∈# mset (prime_factorization_nat n)" for p
        using pf(1) that by simp
      let ?l = "i∈#prime_factorization n. i"
      let ?r = "i∈#mset (prime_factorization_nat n). i"
      show "prod_mset (mset (prime_factorization_nat n)) = normalize n"
        by (simp add: pf(2) prod_mset_prod_list)
    qed
  qed
qed

lemma multiset_prime_factorization_code[code_unfold]: 
  "prime_factorization = (λn. mset (prime_factorization_nat n))"
  by (intro ext multiset_prime_factorization_nat_correct)

lemma divisors_nat: 
  "n  0  set (divisors_nat n) = {p. p dvd n}" "distinct (divisors_nat n)" "divisors_nat 0 = []"
proof -
  show "distinct (divisors_nat n)" "divisors_nat 0 = []" unfolding divisors_nat_def by auto
  assume n: "n  0"
  from n have "n > 0" by auto
  {
    fix x    
    have "(x dvd n) = (x  0  (p. multiplicity p x  multiplicity p n))"
    proof (cases "x = 0")
      case False
      with n > 0 show ?thesis by (auto simp: dvd_multiplicity_eq)
    next
      case True
      with n show ?thesis by auto
    qed
  } note dvd = this
  let ?dn = "set (divisors_nat n)"
  let ?mf = "λ (n :: nat). prime_factorization n"
  have "?dn = prod_list ` set (subseqs (prime_factorization_nat n))" unfolding divisors_nat_def
    using n by auto
  also have " = prod_mset ` mset ` set (subseqs (prime_factorization_nat n))"
    by (force simp: prod_mset_prod_list)
  also have "mset ` set (subseqs (prime_factorization_nat n))
    = { ps. ps ⊆# mset (prime_factorization_nat n)}" 
    unfolding multiset_of_subseqs by simp
  also have " = { ps. ps ⊆# ?mf n}"
    thm multiset_prime_factorization_code[symmetric]
    unfolding multiset_prime_factorization_nat_correct[symmetric] by auto
  also have "prod_mset `  = {p. p dvd n}" (is "?l = ?r")
  proof -
    {
      fix x
      assume "x dvd n"
      from this[unfolded dvd] have x: "x  0" by auto
      from x dvd n x  0 n  0 have sub: "?mf x ⊆# ?mf n"
        by (subst prime_factorization_subset_iff_dvd) auto
      have "prod_mset (?mf x) = x" using x
        by (simp add: prime_factorization_nat)
      hence "x  ?l" using sub by force
    } 
    moreover
    {
      fix x
      assume "x  ?l"
      then obtain ps where x: "x = prod_mset ps" and sub: "ps ⊆# ?mf n" by auto
      have "x dvd n" using prod_mset_subset_imp_dvd[OF sub] n x by simp
    }
    ultimately show ?thesis by blast
  qed
  finally show "set (divisors_nat n) = {p. p dvd n}" .
qed

lemma divisors_int_pos: "x  0  set (divisors_int_pos x) = {i. i dvd x  i > 0}" "distinct (divisors_int_pos x)"
  "divisors_int_pos 0 = []"
proof -
  show "divisors_int_pos 0 = []" by code_simp
  show "distinct (divisors_int_pos x)" 
    unfolding divisors_int_pos_def using divisors_nat(2)[of "nat (abs x)"]
    by (simp add: distinct_map inj_on_def)
  assume x: "x  0"
  let ?x = "nat (abs x)"
  from x have xx: "?x  0" by auto
  from x have 0: " y. y dvd x  y  0" by auto
  have id: "int ` {p. int p dvd x} = {i. i dvd x  0 < i}" (is "?l = ?r")
  proof -
    {
      fix y
      assume "y  ?l"
      then obtain p where y: "y = int p" and dvd: "int p dvd x" by auto
      have "y  ?r" unfolding y using dvd 0[OF dvd] by auto
    }
    moreover
    {
      fix y
      assume "y  ?r"
      hence dvd: "y dvd x" and y0: "y > 0" by auto
      define n where "n = nat y"
      from y0 have y: "y = int n" unfolding n_def by auto
      with dvd have "y  ?l" by auto
    }
    ultimately show ?thesis by blast
  qed
  from xx show "set (divisors_int_pos x) = {i. i dvd x  i > 0}"
    by (simp add: divisors_int_pos_def divisors_nat id) 
qed

lemma divisors_int: "x  0  set (divisors_int x) = {i. i dvd x}" "distinct (divisors_int x)"
  "divisors_int 0 = []"
proof -
  show "divisors_int 0 = []" by code_simp
  show "distinct (divisors_int x)"
  proof (cases "x = 0")
    case True 
    show ?thesis unfolding True by code_simp
  next
    case False
    from divisors_int_pos(1)[OF False] divisors_int_pos(2)
    show ?thesis unfolding divisors_int_def Let_def distinct_append distinct_map inj_on_def by auto
  qed
  assume x: "x  0"
  show "set (divisors_int x) = {i. i dvd x}"
  unfolding divisors_int_def Let_def set_append set_map divisors_int_pos(1)[OF x] using x
  by auto (metis (no_types, lifting) dvd_mult_div_cancel image_eqI linorder_neqE_linordered_idom 
  mem_Collect_eq minus_dvd_iff minus_minus mult_zero_left neg_less_0_iff_less)
qed

definition divisors_fun :: "('a  ('a :: {comm_monoid_mult,zero}) list)  bool" where
  "divisors_fun df  ( x. x  0  set (df x) = { d. d dvd x })  ( x. distinct (df x))"

lemma divisors_funD: "divisors_fun df  x  0  d dvd x  d  set (df x)" 
  unfolding divisors_fun_def by auto

definition divisors_pos_fun :: "('a  ('a :: {comm_monoid_mult,zero,ord}) list)  bool" where
  "divisors_pos_fun df  ( x. x  0  set (df x) = { d. d dvd x  d > 0})  ( x. distinct (df x))"

lemma divisors_pos_funD: "divisors_pos_fun df  x  0  d dvd x  d > 0  d  set (df x)" 
  unfolding divisors_pos_fun_def by auto

lemma divisors_fun_nat: "divisors_fun divisors_nat"
  unfolding divisors_fun_def using divisors_nat by auto

lemma divisors_fun_int: "divisors_fun divisors_int"
  unfolding divisors_fun_def using divisors_int by auto

lemma divisors_pos_fun_int: "divisors_pos_fun divisors_int_pos"
  unfolding divisors_pos_fun_def using divisors_int_pos by auto

end