Theory Prime_Factorization

theory Prime_Factorization
imports Primes Missing_Multiset
(*  
    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 Divides.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 Divides.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: "Divides.divmod_nat n i = (n',mo)" by force
  hence n': "n' = n div i" and mo: "mo = n mod i" by (auto simp: 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.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 "(∏p∈insert 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: "Divides.divmod_nat n i' = (n',m)" by force
    hence n': "n' = n div i'" and m: "m = n mod i'" by (auto simp: 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'':
  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