Theory Wieferich_Kempner

(*  Author:      Jamie Chen
    License:     LGPL      *)

theory "Wieferich_Kempner"
  imports  
  "HOL-Number_Theory.Cong"
  "HOL.Real"
  "HOL.NthRoot"
  "HOL.Transcendental"
  "HOL-Library.Code_Target_Nat"
  "Three_Squares.Three_Squares"
  Main
begin

fun sumpow :: "nat ⇒ nat list ⇒ nat" where [code_computation_unfold, coercion_enabled]:
  "sumpow p l = fold (+) (map (λx. x^p) l) 0"
declare sumpow.simps[code]

definition is_sumpow :: "nat ⇒ nat ⇒ nat ⇒ bool"
  where "is_sumpow p n m ≡ ∃ l. length l = n ∧ m = sumpow p l"

section ‹Technical Lemmas›
text‹We show four lemmas used in the main theorem.›
subsection ‹ Lemma 2.1 in \cite{nathanson1996}›

lemma sum_of_6_cubes:
  fixes A m :: nat
  assumes mLessASq: "m ≤ A2"
  assumes mIsSum3Sq: "is_sumpow 2 3 m"
  shows "is_sumpow 3 6 (6 * A * (A2 + m))"
proof -
  from mIsSum3Sq obtain l where "length l = 3" and " m = sumpow 2 l"
    using is_sumpow_def by blast
  then obtain m1 m2 m3 where "l = [m1, m2, m3]"
    by (smt (verit, best) Suc_length_conv length_0_conv numeral_3_eq_3)

  obtain il where "il = map (int) l" by auto
  obtain iA where "iA = (int A)" by auto
  obtain im where "im = (int m)" by auto

  have "fold (+) (map power2 l) 0 = m12 + m22 + m32"
    using ‹l = [m1, m2, m3]› by simp
  hence "m = m12 + m22 + m32"
    using ‹length l = 3› ‹m = sumpow 2 l› sumpow.simps by presburger
  obtain im1 im2 im3 where "il = [im1, im2, im3]"
    by (simp add: ‹il = map int l› ‹l = [m1, m2, m3]›)
  hence "im1 = int m1" and "im2 = int m2" and "im3 = int m3"
    using ‹il = map int l› ‹l = [m1, m2, m3]› by auto

  obtain x1 x2 x3 where "[x1,x2,x3] = map (λx. A+x) l"
    by (simp add: ‹l =[m1,m2,m3]›)
  obtain x4 x5 x6 where "[x4,x5,x6] = map (λx. A-x) l"
    by (simp add: ‹l =[m1,m2,m3]›)
  obtain ns where "ns =  [x1, x2, x3, x4, x5, x6]" by auto

  have "∀ x ∈ set l. x ≥ 0"
    by auto
  hence "∀ x ∈ set l. x2 ≤ m"
    using‹l = [m1, m2, m3]› ‹m = m12 + m22 + m32› power2_nat_le_imp_le by auto
  have "∀ x ∈ set il. x ≥ 0"
    by (simp add: ‹il = map int l›)
  hence "∀ x ∈ set il. x2 ≤ im"
    using‹im = int m› ‹il = map int l›‹l = [m1, m2, m3]› ‹m = m12 + m22 + m32› power2_nat_le_imp_le 
    by auto

  have "∀ x ∈ set il. iA + x ≥ 0"
    using ‹iA = int A› ‹il = map int l› by auto
  have "int (x1) = iA + im1"
    using ‹[x1, x2, x3] = map ((+) A) l› ‹iA = int A› ‹il = [im1, im2, im3]› ‹il = map int l› by force
  have "int (x2) = iA + im2"
    using ‹[x1, x2, x3] = map ((+) A) l› ‹iA = int A› ‹il = [im1, im2, im3]› ‹il = map int l› by force
  have "int (x3) = iA + im3"
    using ‹[x1, x2, x3] = map ((+) A) l› ‹iA = int A› ‹il = [im1, im2, im3]› ‹il = map int l› by force
  have "A ≥ m1"
    by (metis ‹m = m12 + m22 + m32› add_leE mLessASq power2_nat_le_eq_le)
  hence "iA - im1 ≥ 0"
    using ‹iA = int A› ‹il = [im1, im2, im3]› ‹il = map int l› ‹l = [m1, m2, m3]› by auto
  hence "int x4 = iA - im1"
    using ‹[x4, x5, x6] = map ((-) A) l› ‹iA = int A› ‹im1 = int m1› ‹l = [m1, m2, m3]› by auto
  have "A ≥ m2"
    by (metis ‹m = m12 + m22 + m32› add_leE mLessASq power2_nat_le_eq_le)
  hence "iA - im2 ≥ 0"
    using ‹iA = int A› ‹il = [im1, im2, im3]› ‹il = map int l› ‹l = [m1, m2, m3]› by auto
  hence "int x5 = iA - im2"
    using ‹[x4, x5, x6] = map ((-) A) l› ‹iA = int A› ‹im2 = int m2›‹l = [m1, m2, m3]› by auto
  have "A ≥ m3"
    by (metis ‹m = m12 + m22 + m32› add_leE mLessASq power2_nat_le_eq_le)
  hence "iA - im3 ≥ 0"
    using ‹iA = int A› ‹il = [im1, im2, im3]› ‹il = map int l› ‹l = [m1, m2, m3]› by auto
  hence "int x6 = iA - im3"
    using ‹[x4, x5, x6] = map ((-) A) l› ‹iA = int A› ‹im3 = int m3› ‹l = [m1, m2, m3]› by auto

  have "6*A*(A2+m) = 6*A*(A2 + sumpow 2 l)"
    by (simp add:‹m = sumpow 2 l›)
  have distr: "… = 6*A^3 +6*A*(sumpow 2 l)"
    by (simp add: distrib_left power2_eq_square power3_eq_cube) 
  have "… = 6*A*(A2 + m12 + m22 + m32)"
    using ‹m = m12 + m22 + m32› ‹m = sumpow 2 l› by (metis distr group_cancel.add1)
  hence expanded: "int … = 6*iA*(iA2 +  im12 +im22 + im32)"
    using ‹iA = int A› ‹im1 = int m1› ‹im2 = int m2› ‹im3 = int m3› by simp
  have sixcubes: "… = (iA + im1)^3 + (iA - im1)^3 +
                       (iA + im2)^3 + (iA - im2)^3 +
                       (iA + im3)^3 + (iA - im3)^3"
    by Groebner_Basis.algebra
  have "… = int x1^3 + int x2^3 + int x3^3 + int x4^3 + int x5^3 + int x6^3"
    using ‹int x1 = iA + im1› ‹int x2 = iA + im2› ‹int x3 = iA + im3› 
          ‹int x4 = iA - im1› ‹int x5 = iA - im2› ‹int x6 = iA - im3› by simp
  hence "6*A*(A2+m) = x1^3 + x2^3 + x3^3 + x4^3 + x5^3 + x6^3"
    using distr ‹m = sumpow 2 l› ‹m = m12 + m22 + m32› expanded of_nat_eq_iff sixcubes 
    by (smt (verit) of_nat_add of_nat_power)
  have "map (λx. x^3) ns = [x1^3, x2^3, x3^3, x4^3, x5^3, x6^3]"
    by (simp add: ‹ns = [x1, x2, x3, x4, x5, x6]›)
  hence "sumpow 3 ns = x1^3 + x2^3 + x3^3 + x4^3 + x5^3 + x6^3"
    by (simp add: ‹ns = [x1, x2, x3, x4, x5, x6]›)
  hence "6*A*(A2+m) = sumpow 3 ns"
    using ‹6 * A * (A2 + m) = x1 ^ 3 + x2 ^ 3 + x3 ^ 3 + x4 ^ 3 + x5 ^ 3 + x6 ^ 3› by presburger
  hence "length ns = 6 ∧ 6*A*(A2+m) = sumpow 3 ns"
    by (simp add: ‹ns = [x1, x2, x3, x4, x5, x6]›)
  then show ?thesis
    using is_sumpow_def by blast
qed

subsection ‹ Lemma 2.2 in \cite{nathanson1996}›

lemma if_cube_cong_then_cong:
  fixes t :: nat
  fixes b1 b2 :: int
  assumes odd: "odd b1 ∧ odd b2"
  assumes "b1 > 0 ∧ b2 > 0"
  assumes tGeq1: "t ≥ 1"
  assumes "[b1^3 = b2^3] (mod 2^t)"
  shows "[b1 ≠ b2] (mod 2^t) ⟹ [b1^3 ≠ b2^3] (mod 2^t)"
proof -
  have "[b2^3 - b1^3 = 0] (mod 2^t)"
    using ‹[b1 ^ 3 = b2 ^ 3] (mod 2 ^ t)› cong_diff_iff_cong_0 cong_sym_eq by blast
  have "(b2 - b1)*(b22 + b2*b1 + b12) = b2^3 - b1^3"
    by Groebner_Basis.algebra
  hence "[(b2-b1)*(b22 + b2*b1 + b12) = 0] (mod 2^t)"
    using ‹[b2 ^ 3 - b1 ^ 3 = 0] (mod 2 ^ t)› by auto
  have "coprime (b22 + b2*b1 + b12) (2^t)"
    by (simp add: odd)
  hence "[b2 - b1 = 0] (mod 2^t)"
    by (metis ‹[(b2 - b1) * (b22 + b2 * b1 + b12) = 0] (mod 2 ^ t)› cong_mult_rcancel mult_zero_left)
  hence "[b1 = b2] (mod 2^t)"
    using cong_diff_iff_cong_0 cong_sym_eq by blast
  assume "[b1 ≠ b2] (mod 2^t)"
  then show ?thesis
    using ‹[b1 = b2] (mod 2^t)› by auto
qed

lemma every_odd_nat_cong_cube:
  fixes t w :: nat
  assumes tPositive: "t ≥ 1"
  assumes wOdd: "odd w"
  shows "∃b. odd b ∧ [w = b^3] (mod 2^t)"
proof (rule ccontr)
  assume "¬?thesis"
  hence "∀b. odd b ⟶ [w ≠ b^3] (mod 2^t)"
    by blast
  obtain b::nat where "odd b"
    using odd_numeral by blast
  hence "[w ≠ b^3] (mod 2^t)"
    by (simp add: ‹∀b. odd b ⟶ [w ≠ b ^ 3] (mod 2 ^ t)›)
  obtain bSet::"nat set" where bSet_def: "bSet = {x. odd x ∧ x < 2^t}"
    by auto
  obtain w' where w'_def: "w' = w mod 2^t"
    by auto
  hence "odd w'"
    by (metis dvd_mod_imp_dvd dvd_power order_less_le_trans tPositive wOdd zero_less_one)
  have "w' < 2^t"
    by (simp add: ‹w' = w mod 2^t›)
  hence "w' ∈ bSet"
    by (simp add: ‹bSet = {x. odd x ∧ x < 2 ^ t}› ‹odd w'›)
  define bSetMinusW'::"nat set" where "bSetMinusW' = {x. x ∈ bSet ∧ x ≠ w'}"

  have "∀x ∈ bSet. [w ≠ x^3] (mod 2^t)"
    using ‹∀b. odd b ⟶ [w ≠ b ^ 3] (mod 2 ^ t)› ‹bSet = {x. odd x ∧ x < 2 ^ t}› by blast
  have "∀x ∈ bSet. odd (x^3)"
    using bSet_def by simp
  have "even (2^t)"
    using tPositive by fastforce
  hence "∀x ∈ bSet. odd (x^3 mod 2^t)"
    using ‹∀x∈bSet. odd (x^3)› dvd_mod_iff by blast
  hence "∀ x ∈ bSet. (x^3 mod 2^t) ∈ bSet"
    by (simp add: ‹bSet = {x. odd x ∧ x < 2 ^ t}›)
  hence "∀ x ∈ bSet. (x^3 mod 2^t) ∈ bSetMinusW'"
    using ‹∀x∈bSet. [w ≠ x ^ 3] (mod 2 ^ t)› bSetMinusW'_def w'_def cong_def
    by (metis (mono_tags, lifting) mem_Collect_eq unique_euclidean_semiring_class.cong_def)
  hence "∀ x ∈ bSet. ∃ y ∈ bSetMinusW'. [y = x^3] (mod 2^t)"
    using cong_mod_right cong_refl by blast
  hence "∀ x ∈ bSet. ∃! y ∈ bSetMinusW'. [y = x^3] (mod 2^t)"
    by (metis (no_types, lifting) bSet_def bSetMinusW'_def unique_euclidean_semiring_class.cong_def 
        mem_Collect_eq mod_less)
  have "∀ x1 ∈ bSet.∀ x2 ∈ bSet. odd x1 ∧ odd x2"
    using bSet_def by auto

  hence "∀ x1 ∈ bSet. ∀ x2 ∈ bSet. [x1 ≠ x2] (mod 2^t) ⟶ [x1^3 ≠ x2^3] (mod 2^t)"
    by (metis (mono_tags) cong_int_iff even_of_nat if_cube_cong_then_cong odd_pos of_nat_0_less_iff 
        of_nat_numeral of_nat_power tPositive)
  hence "∀ x1 ∈ bSet. ∀ x2 ∈ bSet. x1 ≠ x2 ⟶ [x1^3 ≠ x2^3] (mod 2^t)"
    using bSet_def cong_less_modulus_unique_nat by blast
  hence "inj_on (λ x. x^3 mod 2^t) bSet"
    by (meson inj_onI unique_euclidean_semiring_class.cong_def)

  have "bSetMinusW' ⊂ bSet"
    using bSetMinusW'_def ‹w' ∈ bSet› by blast 
  moreover hence "card bSetMinusW' < card bSet"
    by (simp add: bSet_def psubset_card_mono)
  moreover have "(λ x. x^3 mod 2^t) ` bSet ⊆ bSetMinusW'"
    using ‹∀x∈bSet. x ^ 3 mod 2 ^ t ∈ bSetMinusW'› by blast
  ultimately have "card ((λ x. x^3 mod 2^t) ` bSet) < card bSet"
    by (metis card.infinite less_zeroE psubset_card_mono subset_psubset_trans)
  hence "¬inj_on (λ x. x^3 mod 2^t) bSet"
    using pigeonhole by auto
  thus False
    using ‹inj_on (λx. x ^ 3 mod 2 ^ t) bSet› by blast
qed

subsection ‹ Lemma 2.3 in \cite{nathanson1996}›
text‹It is this section in which we use the Three Squares Theorem AFP Entry \cite{Three_Squares-AFP}.›

lemma sum_of_3_squares_exceptions:
  fixes m::nat
  assumes notSum3Sq: "¬is_sumpow 2 3 m"
  shows "6*m mod 96 ∈ {0,72,42,90}"
proof -
  have "¬(∃ l. length l = 3 ∧ m = sumpow 2 l)"
    using is_sumpow_def notSum3Sq by blast
  hence "¬(∃ a b c. m = sumpow 2 [a,b,c])"
    by (metis length_Cons list.size(3) numeral_3_eq_3)
  hence "¬(∃ a b c. m = fold (+) (map power2 [a,b,c]) 0)"
    by auto
  hence "¬(∃ a b c. m = fold (+) [a2,b2,c2] 0)"
    by auto
  hence "¬(∃ a b c. m = a2 + b2 + c2)"
    by (simp add: add.commute)
  then obtain s t where "(m = 4^s*(8*t + 7))"
    using three_squares_iff by auto
  obtain h'::"nat set" where "h' = {0,72,42,90}" 
    by auto
  consider (s0) "s = 0"| (s1) "s = 1" | (s2) "s≥2"
    by arith
  hence "6*4^s*(8*t + 7) mod 96 ∈ h'"
  proof cases
    case s0
    consider (even) "even t" | (odd) "odd t"
      by auto
    thus "6*4^s*(8*t + 7) mod 96 ∈ h'"
    proof cases
      case even
      obtain t' where "t' =  t div 2"
        by simp
      hence "6*4^s*(8*t+7) = 96*t' + 42"
        using s0 even by fastforce
      hence "6*4^s*(8*t+7) mod 96 = 42"
        by (simp add: cong_0_iff cong_add_rcancel_0_nat)
      have "42 ∈ h'"
        by (simp add: ‹h' = {0, 72, 42, 90}›)
      thus "6*4^s*(8*t + 7) mod 96 ∈ h'"
        by (simp add: ‹6 * 4 ^ s * (8 * t + 7) mod 96 = 42›)
    next
      case odd
      obtain t' where "t' =  (t-1) div 2"
        by simp
      hence "6*4^s*(8*t+7) = 6*(8*(2*t'+1)+7)"
        by (simp add: s0 odd)
      hence "6*4^s*(8*t+7) mod 96 = 90"
        using cong_le_nat by auto
      have "90 ∈ h'"
        by (simp add: ‹h' = {0, 72, 42, 90}›)
      thus "6*4^s*(8*t + 7) mod 96 ∈ h'"
        by (simp add: ‹6 * 4 ^ s * (8 * t + 7) mod 96 = 90›)
    qed
  next
    case s1
    have "6*4^s*(8*t+7) = 96*(2*t + 1) + 72"
      using s1 by auto
    hence "[6*4^s*(8*t+7) = 72] (mod 96)"
    by (metis cong_add_rcancel_0_nat cong_mult_self_left)
  hence "6*4^s*(8*t+7) mod 96 = 72 mod 96"
    using unique_euclidean_semiring_class.cong_def by blast
    hence "6*4^s*(8*t+7) mod 96 = 72"
      by auto
    have "72 ∈ h'"
      by (simp add: ‹h' = {0, 72, 42, 90}›)
    thus "6*4^s*(8*t + 7) mod 96 ∈ h'"
      by (simp add:‹6 * 4 ^ s * (8 * t + 7) mod 96 = 72›)
  next
    case s2
    obtain s' where "s' = s-2"
      by simp
    have "6*4^s*(8*t+7) = 6*4^2*4^s'*(8*t+7)"
      by (metis ‹s' = s - 2› le_add_diff_inverse mult.assoc power_add s2)
    have "… = 96*4^s'*(8*t+7)"
      by auto
    hence "[6*4^s*(8*t+7) = 0] (mod 96)"
      by (metis ‹6*4^s*(8*t + 7) = 6*42*4^s'*(8*t + 7)› cong_modulus_mult_nat cong_mult_self_left)
    hence "6*4^s*(8*t+7) mod 96 = 0 mod 96"
    using unique_euclidean_semiring_class.cong_def by blast
    hence "6*4^s*(8*t+7) mod 96 = 0"
      by auto
    have "0 ∈ h'"
      by (simp add: ‹h' = {0, 72, 42, 90}›)
    thus "6*4^s*(8*t + 7) mod 96 ∈ h'"
      by (simp add: ‹6 * 4 ^ s * (8 * t + 7) mod 96 = 0›)
  qed
  thus ?thesis
     by (metis ‹h' = {0, 72, 42, 90}› ‹m = 4 ^ s * (8 * t + 7)› mult.assoc)
qed

lemma values_geq_22_cubed_can_be_normalised:
  fixes r :: nat
  assumes rLarge: " r ≥ 10648" (*10648 = 223*)
  obtains d m where "d ≥ 0" and "d ≤ 22" and "r = d^3 + 6*m" and "is_sumpow 2 3 m"
proof -
  define MultiplesOf6LessThan96::"nat set" where 
    "MultiplesOf6LessThan96 = {0,6,12,18,24,30,36,42,48,54,60,66,72,78,84,90}"
  have "∀m' ∈  MultiplesOf6LessThan96. m' < 96"
    unfolding MultiplesOf6LessThan96_def by auto
  have "∀ m' ∈ MultiplesOf6LessThan96. m' mod 6 = 0"
    unfolding MultiplesOf6LessThan96_def by auto
  hence "∀m ∈ MultiplesOf6LessThan96. ∃n. m = 6*n"
     by fastforce
  have "∀m' ∈ MultiplesOf6LessThan96. m' mod 96 = m'"
    by (simp add: ‹∀m'∈MultiplesOf6LessThan96. m' < 96›)

  define H'::"nat set" where "H' = {0,42,72,90}"

  have "∀m'. 6*m' mod 96 ∉ H' ⟶ is_sumpow 2 3 m'"
    unfolding H'_def using sum_of_3_squares_exceptions by blast
  hence "∀m'. (6*m' mod 96 ∈ MultiplesOf6LessThan96 - H') ⟶ is_sumpow 2 3 m'"
    by fastforce

  define H::"nat set" where "H = {6, 12, 18, 24, 30, 36, 48, 54, 60, 66, 78, 84}"
  have "H =  MultiplesOf6LessThan96 - H'"
    unfolding H_def MultiplesOf6LessThan96_def H'_def by auto
  hence "∀m'. 6*m' mod 96 ∈ H ⟶ is_sumpow 2 3 m'"
    using ‹∀m'. 6 * m' mod 96 ∈ MultiplesOf6LessThan96 - H' ⟶ is_sumpow 2 3 m'› by blast

  define D::"nat set" where "D = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22}"
  have "∀d ∈ D. d ≥ 0 ∧ d ≤ 22"
    using D_def by auto
  define residue::"nat set" where "residue = {x. ∃ h ∈ H. ∃ d ∈ D. x = (d^3 + h) mod 96}"
(*for some reason Isabelle doesn't like it when you do this in steps larger than 10*)
  have "∀x<10. x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 10. x < 19 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 19. x < 28 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 28. x < 37 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 37. x < 46 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 46. x < 55 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 55. x < 64 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 64. x < 73 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 73. x < 82 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 82. x < 91 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  moreover have "∀x≥ 91. x < 96 ⟶ x ∈ residue"
    unfolding residue_def H_def D_def by auto
  ultimately have "∀x < 96. x ∈ residue"
    by (metis leI)

  have "∀ r ∈ residue. ∃ h ∈ H. ∃ d ∈ D. r = (d^3 + h) mod 96"
    unfolding residue_def by auto
  have "∀ h ∈ H. ∀ d ∈ D. (d^3 + h) mod 96 ∈ residue"
    by (simp add: ‹∀x<96. x ∈ residue›)
  have "r mod 96 ∈ residue"
    by (simp add: ‹∀x<96. x ∈ residue›)
  have "∀d ∈ D. r ≥ d^3"
    unfolding D_def using rLarge by auto
  hence "∃ h ∈ H. ∃ d ∈ D. [d^3 + h = r] (mod 96)"
    using ‹r mod 96 ∈ residue›
    by (metis ‹∀r∈residue. ∃h∈H. ∃d∈D. r = (d ^ 3 + h) mod 96› unique_euclidean_semiring_class.cong_def)
  hence "∃ h ∈ H. ∃ d ∈ D. [r - d^3 = h] (mod 96)"
    by (smt (z3) ‹∀d∈D. d^3 ≤ r› add.commute add_diff_cancel_right' cong_diff_nat cong_refl le_add2)
  then obtain h d where "h ∈ H" and "d ∈ D" and "[r - d^3 = h] (mod 96)"
    by auto
  then obtain h' where "h' mod 96 = h" and "r - d^3 = h'"
    using ‹H = MultiplesOf6LessThan96 - H'› ‹∀m'∈MultiplesOf6LessThan96. m' mod 96 = m'› 
    by (simp add: unique_euclidean_semiring_class.cong_def)
  have "[h = 0] (mod 6)"
    using ‹H = MultiplesOf6LessThan96 - H'› ‹∀m∈MultiplesOf6LessThan96. ∃n. m = 6*n› ‹h ∈ H› 
      cong_mult_self_left by fastforce
  hence "[h' = 0] (mod 6)"
    using ‹[r-d^3 = h] (mod 96)› ‹r-d^3 = h'› unique_euclidean_semiring_class.cong_def
    by (metis cong_dvd_modulus_nat dvd_triv_right num_double numeral_mult)
  then obtain m where "6*m = h'"
    by (metis cong_0_iff dvd_def)
  moreover hence "is_sumpow 2 3 m"
    using ‹∀m'. 6 * m' mod 96 ∈ H ⟶ is_sumpow 2 3 m'› by (simp add:‹h ∈ H› ‹h' mod 96 = h›)
  moreover have "d ≥ 0" 
    by auto
  moreover have "d ≤ 22"
    by (simp add: ‹∀d∈D. 0 ≤ d ∧ d ≤ 22› ‹d ∈ D›)
  ultimately show "?thesis"
    by (metis ‹∀d∈D. d ^ 3 ≤ r› ‹d ∈ D›‹r - d ^ 3 = h'› le_add_diff_inverse that)
qed

subsection ‹Lemma 2.4 in \cite{nathanson1996}›

partial_function(tailrec) list_builder :: "nat ⇒ nat ⇒ nat list ⇒ nat list" 
  where [code_computation_unfold, coercion_enabled]:
    "list_builder m n l = (if n = 0 then l else (list_builder m (n-1) (m#l)))"
declare list_builder.simps[code]

partial_function(tailrec) dec_list :: "nat ⇒ nat list ⇒ nat list" 
  where [code_computation_unfold, coercion_enabled]:
    "dec_list depth l = (if (tl l) = [] then list_builder (hd l + 1) (depth+1) [] else 
    if hd l ≤ hd (tl l) + 1 then dec_list (depth+1) ((hd (tl l) + 1)#(tl (tl l))) else
    list_builder (hd (tl l) + 1) (depth + 2) (tl (tl l)))"
declare dec_list.simps[code]

partial_function(tailrec) sumcubepow_finder ::"nat ⇒ nat list ⇒ nat list" 
  where [code_computation_unfold, coercion_enabled]:
    "sumcubepow_finder n l = (if (sumpow 3 l < n) then 
    (sumcubepow_finder n ((Suc (hd l))#(tl l)))
    else if (sumpow 3 l) = n then l else sumcubepow_finder n  (dec_list 0 l))"
declare sumcubepow_finder.simps[code]

lemma leq_40000_is_sum_of_9_cubes:
  fixes n :: nat
  assumes "n ≤ 40000"
  shows "is_sumpow 3 9 n" and "n > 8042 ⟶ is_sumpow 3 6 n"
proof -
  (* split computation into chunks so that computation is done in parallel *)
  let ?A1 = "[8043 ..< 15000]" 
  let ?A2 = "[15000 ..< 23000]" 
  let ?A3 = "[23000 ..< 29000]" 
  let ?A4 = "[29000 ..< 35000]" 
  let ?A5 = "[35000 ..< 40001]" 
  (* let expression to avoid duplicate computation of result *)
  let ?test = "λ n. let result = sumcubepow_finder n [0,0,0,0,0,0] in
     sumpow 3 result = n ∧ length result = 6" 
  have split: "⋀ n. n ≤ 40000 ⟹ n > 8042 ⟹ 
    n ∈ set ?A1 ∨ n ∈ set ?A2 ∨ n ∈ set ?A3 ∨ n ∈ set ?A4 ∨ n ∈ set ?A5" 
    unfolding set_upt by simp linarith
  have "∀ n ∈ set ?A1. ?test n"
    by eval (*takes a while to compute*)
  moreover have "∀n ∈ set ?A2. ?test n"
    by eval (*takes a while to compute*)
  moreover have "∀n ∈ set ?A3. ?test n"
    by eval (*takes a while to compute*)
  moreover have "∀n ∈ set ?A4. ?test n"
    by eval (*takes a while to compute*)
  moreover have "∀n ∈ set ?A5. ?test n"
    by eval (*takes a while to compute*)
  ultimately have "∀n ≤ 40000. n > 8042 ⟶ ?test n" using split by blast
  then have "∀n≤40000. n > 8042 ⟶ is_sumpow 3 6 n"
    using is_sumpow_def by metis
  then show "n > 8042 ⟶ is_sumpow 3 6 n"
    using ‹n ≤ 40000› by simp
  have "∀n≤40000. n > 8042 ⟶ (∃ l. length l = 6 ∧ n = sumpow 3 l)"
    using is_sumpow_def ‹∀n≤40000. n > 8042 ⟶ is_sumpow 3 6 n› by simp
  then have "∀n≤40000. n > 8042 ⟶ (∃ l. length (l@[0,0,0]) = 9 ∧ n = sumpow 3 (l@[0,0,0]))"
    by simp
  then have "∀n≤40000. n > 8042 ⟶ is_sumpow 3 9 n"
    using is_sumpow_def by blast
  have "∀n ≤ 8042. let res = sumcubepow_finder n [0,0,0,0,0,0,0,0,0] in 
      sumpow 3 res = n ∧ length res = 9"
    by eval
  then have "∀n≤8042. is_sumpow 3 9 n"
    using is_sumpow_def by metis
  then have "∀n ≤ 40000. is_sumpow 3 9 n"
    using ‹∀n≤40000. n > 8042 ⟶ is_sumpow 3 9 n› by auto
  then show "is_sumpow 3 9 n"
    using ‹ n ≤ 40000› by simp
qed

section ‹Wieferich--Kempner Theorem›
text‹Theorem 2.1 in \cite{nathanson1996}›

theorem Wieferich_Kempner:
  fixes N :: nat
  shows "is_sumpow 3 9 N"
proof cases
  consider (ge8pow10) "N > 8^10" | (leq8pow10) "N ≤ 8^10"
    by arith
  thus ?thesis
  proof cases
    case ge8pow10
    define n::int where "n = ⌊root 3 N⌋"
    hence "n ≥ ⌊root 3 (8^10)⌋"
      by (meson floor_mono ge8pow10 less_imp_le_nat numeral_power_le_of_nat_cancel_iff 
          real_root_le_iff zero_less_numeral)
    have "(2::nat)^3 = 8"
      by simp
    hence "8^10 = ((2::nat)^10)^3"
      by simp
    have "root 3 ((2::nat)^10)^3 = (2::nat)^10"
      by simp
    hence "n ≥ ⌊(2::nat)^10⌋"
      by (metis ‹8^10 = (2^10)^3› ‹⌊root 3 (8^10)⌋ ≤ n› of_nat_numeral of_nat_power power3_eq_cube 
          real_root_mult)
    hence ngeq2pow10: " n ≥ (2::nat)^10"
      by auto
    obtain k::int where "k = ⌈(log 8 (N/8))/3⌉-1"
      by simp
    hence "3*k < (log 8 (N/8))"
      by linarith
    hence "8*(8 powr (3*k)) < N"
      using ge8pow10 less_log_iff by simp
    have "k+1 ≥ (log 8 (N/8))/3"
      by (simp add: ‹k = ⌈log 8 (real N / 8) / 3⌉ - 1›)
    hence "8 powr (3*(k+1)) ≥ 8 powr (log 8 (N/8))"
      by simp
    hence "8*(8 powr (3*(k+1))) ≥ N"
      using ge8pow10 by simp
    have "(N/8) > 8^9"
      using ge8pow10 by simp
    hence "(log 8 (N/8)) > 9"
      by (metis less_log_of_power numeral_One numeral_less_iff of_nat_numeral semiring_norm(76))
    hence "(log 8 (N/8))/3 > 3"
      by simp
    hence "k ≥ 3"
      using ‹k = ⌈(log 8 (N/8))/3⌉-1› by simp

    define i::int where "i = ⌊root 3 (N - (8*(8 powr (3*k))))⌋"
    hence "i ≤ root 3 (N - (8*(8 powr (3*k))))"
      by fastforce
    have "root 3 (N - (8*(8 powr (3*k)))) ≥ 0"
      using ‹8*(8 powr (3 * k)) < N› by force
    hence "i ≥ 0"
      by (simp add: i_def)
    hence "i^3 ≤ (root 3 (N - (8*(8 powr (3*k)))))^3"
      using ‹i ≤ root 3 (N - 8 * 8 powr (3*k))› power_mono by fastforce    
    hence "i^3 ≤ N - (8*(8 powr (3*k)))"
      using ‹8 * 8 powr (3*k) < N› by force
    hence exp01:"8*(8 powr (3*k)) ≤ N - i^3"
      by simp
    have "i+1 > root 3 (N - (8*(8 powr (3*k))))"
      by (simp add: i_def)
    hence "(i+1)^3 > (root 3 (N - (8*(8 powr (3*k)))))^3"
      by (metis ‹0 ≤ root 3 (N - 8 * 8 powr (3*k))› of_int_power power_strict_mono zero_less_numeral)
    hence "(i+1)^3 > (N - (8*(8 powr (3*k))))"
      using ‹8 * 8 powr (3*k) < N› by force
    hence exp02:"8*(8 powr (3*k)) > N - (i+1)^3"
      by simp

    have "i ≥ 1"
    proof (rule ccontr)
      assume "¬ i≥ 1"
      hence "i = 0"
        using ‹i ≥ 0› by simp
      hence "8*(8 powr (3*k)) ≤ N"
        using ‹8*(8 powr (3*k)) ≤ N - i^3› by simp
      moreover have "8*(8 powr (3*k)) + 1 > N"
        using ‹i = 0› ‹(N - (i+1)^3) < 8*(8 powr (3*k))› ge8pow10 by auto
      moreover have "∀ x. x ≤ N ∧ x + 1 > N ⟶ x = N"
        by linarith
      have "8*8^(nat(3*k)) = 8*(8 powr (3*k))"
        by (smt (verit, ccfv_SIG) ‹3 ≤ k› powr_real_of_int)
      ultimately have "8*(8 powr (3*k)) = N"
        by (metis nat_le_real_less nle_le of_nat_le_iff of_nat_numeral of_nat_power power.simps(2))
      thus False
        using ‹8 * 8 powr (3 * k) < N› by simp
    qed
    
    have "8 powr (3*(k+1)) = (8 powr (k+1))^3"
      by (metis of_int_mult of_int_numeral powr_ge_zero powr_numeral powr_powr powr_powr_swap) 
    have "∀(x::int)≥1. x^3 - (x-1)^3  = 3*x^2 - 3*x +1"
      by Groebner_Basis.algebra
    hence "∀(x::int)≥1. x^3 - (x-1)^3  < 3*x^2"
      by auto
    have exp03:"∀(x::int)≥1. x ≤ n ⟶ 3*x^2 ≤ 3*(root 3 N)^2"
      using n_def by (simp add: le_floor_iff)
    have "root 3 N ≤ root 3 (8*(8 powr (3*(k+1))))"
      using ‹N ≤ 8 * 8 powr (3 * (k + 1))› by force
    moreover have "… = (root 3 8)*(root 3 ((8 powr ((k+1)))^3))"
      using ‹8 powr (3*(k+1)) = (8 powr (k+1))^3› real_root_mult by simp
    moreover have "… = 2*(8 powr ((k+1)))"
      using odd_real_root_power_cancel by simp
    ultimately have exp04:"(root 3 N)^2 ≤ (2*(8 powr (k+1)))^2"
      by (metis of_nat_0_le_iff power_mono real_root_pos_pos_le)
    moreover hence exp05: "… = (4*(8 powr ((k+1)))^2)"
      using four_x_squared by presburger
    moreover hence exp06:"… = (4*(8 powr (2*(k+1))))"
      by (smt (verit) of_int_add power2_eq_square powr_add)
    moreover hence "… = (8*(8 powr (2*k + 2))/2)"
      by simp
    moreover hence exp07:"… = (8 powr (2*k + 3)/2)"
      by (simp add: add.commute powr_mult_base)
    ultimately have exp08:"3*(root 3 N)^2 ≤ (3*8 powr (2*k + 3))/2"
      by linarith
    have exp09:"∀ x ≥ 1. x ≤ n ⟶ real_of_int (x^3 - (x-1)^3) < real_of_int (3*x^2)"
      using ‹∀(x::int)≥1. x^3 - (x-1)^3  < 3*x^2› by presburger 
    hence "∀ x ≥ 1. x ≤ n ⟶x^3 - (x-1)^3 < 3*(root 3 N)^2"
      using exp03 less_le_trans 
      by fastforce
    hence exp10:"∀ x ≥ 1. x ≤ n ⟶x^3 - (x-1)^3 ≤ (3*8 powr (2*k + 3))/2"
      using exp08 by fastforce

    have "(root 3 N)^3 = N"
      by simp
    have "root 3 N < n+1"
      using n_def by linarith
    hence "(root 3 N)^3 < (n+1)^3"
      by (metis of_int_power of_nat_0_le_iff real_root_pos_pos_le zero_less_numeral power_strict_mono)
   
    text‹The following few lines have been slightly modified from the original proof to simplify 
formalisation.  This does not affect the proof in any meaningful way.›
    hence exp11:"N - n^3 < (n+1)^3 - n^3"
      using ‹(root 3 N)^3 = N› by linarith
    have exp12:"… = 3*n^2 + 3*n + 1"
      by Groebner_Basis.algebra
    have exp13:"… ≤ 6*n^2 + 1"
      using ‹n ≥ (2::nat)^10› power_strict_mono by (simp add: self_le_power)
    have "… ≤ 6*(root 3 N)^2 + 1"
      using n_def by auto
    have "… ≤ (3*8 powr (2*k + 3))+1"
      using ‹3*(root 3 N)^2 ≤ (3*8 powr (2*k + 3))/2› by auto
    hence exp14:"… ≤ 4*8 powr (2*k + 3)"
      using ‹3 ≤ k› ge_one_powr_ge_zero by auto
    hence exp15:"… ≤ 8*8 powr(3*k)"
      by (smt (verit, best) ‹3 ≤ k› of_int_le_iff powr_mono)

    have "6*n^2≤ 3*8 powr (2*k + 3)"
      using ‹6*(root 3 N)2 + 1 ≤ 3*8 powr (2*k + 3)+1› ‹(6*n2 + 1) ≤ 6*(root 3 N)2 + 1› by linarith
    hence "i ≤ n-1"
      using exp12 exp14 exp13 exp01 i_def n_def exp11 exp15
      by (smt (verit, ccfv_SIG) floor_mono of_int_less_iff real_root_le_iff zero_less_numeral)
    hence "N - (i+1)^3 ≤ 8*8 powr (3*k)"
      using ‹(int N - (i + 1) ^ 3) < 8 * 8 powr (3*k)› by linarith
    hence "N ≥ (i+1)^3"
      using exp04 exp06 exp15 exp07 exp01 exp03 exp09 ‹i ≤ n - 1› ‹0 ≤ i› 
      by (smt (verit, ccfv_threshold) divide_cancel_right exp05 of_int_0_le_iff of_int_diff)
    have exp16:"((i+1)^3 - i^3) ≤ 3*8 powr (2*k + 3)"
      using exp10 exp05 exp04 exp06 exp07 exp03 exp09 ‹i ≤ n-1› ‹0≤i›
      by (smt (verit, ccfv_SIG)  divide_cancel_right powr_non_neg)
    have "(i^3 - (i-1)^3) ≤ 3*8 powr (2*k + 3)"
      using exp10 exp05 exp04 exp06 exp07 exp03 exp09 ‹i ≤ n-1› ‹1≤i›
      by (smt (verit, ccfv_SIG) divide_cancel_right powr_non_neg)
    have "i^3 > (i-1)^3"
      by (smt (verit, best) power_minus_Bit1 power_mono_iff zero_less_numeral)
    have "N - (i-1)^3 = ((N - (i-1)^3)-(N-i^3)) + ((N-i^3) - (N-(i+1)^3)) + (N-(i+1)^3)"
      by simp
    have "… = (i^3 - (i-1)^3) + ((i+1)^3 - i^3) + (N-(i+1)^3)"
      by simp
    have exp17: "… < 3*8 powr (2*k + 3) + 8*8 powr (3*k)"              
      using ‹N - (i+1)^3 ≤ 8*8 powr (3*k)› ‹i ≤ n-1› ‹1≤i› ‹0 ≤ i› exp08 exp10 exp03 exp09
      by (smt (verit) field_sum_of_halves of_int_add of_int_diff power_strict_mono zero_less_numeral)
    have exp18:"… ≤ 11*(8 powr (3*k))"
      by (simp add: ‹3 ≤ k›)
    have "(even (i^3) ∧ odd ((i-1)^3)) ∨ (even ((i-1)^3) ∧ odd (i^3))"
      by simp
    hence "(even (N - i^3) ∧ odd (N - (i-1)^3)) ∨ (even (N - (i-1)^3) ∧ odd (N - i^3))"
      by auto
    hence "(even (N-(i-1)^3)) ⟶ odd (N-i^3)"
      by blast 
    obtain a::nat where a_def: "if odd (N-(i-1)^3) then a = i-1 else a = i"
      by (metis ‹0 ≤ i› ‹1 ≤ i› diff_ge_0_iff_ge nonneg_int_cases)
    have "a = int a"
      by simp
    consider (odd) "odd (N - (i-1)^3)" | (even) "even (N - (i-1)^3)"
      by blast
    hence "odd (N - a^3)"
    proof cases
      case odd
      have "int a = i-1"
        using a_def odd by simp
      have "odd (N - (i-1)^3)"
        using odd by simp
      hence "odd (N - (int a)^3)"
        using ‹a = i-1› by simp
      have "a^3 ≤ N"
        using ‹N ≥ (i+1)^3› ‹int a = i-1› power_mono ‹i ≥ 1›
        by (smt (verit, ccfv_SIG) of_nat_le_of_nat_power_cancel_iff)
      hence "N - a^3 = N - (int a)^3"
        by (simp add: of_nat_diff)
      thus ?thesis
        using ‹odd (N - (int a)^3)› by presburger
    next
      case even
      have "a = i"
        using a_def even by simp
      have "odd (N - i^3)"
        using ‹(even (N-(i-1)^3)) ⟶ odd (N-i^3)› even by simp
      hence "odd (N - (int a)^3)"
        using ‹a = i› by simp
      have "a^3 ≤ N"
        using ‹N ≥ (i+1)^3› ‹int a = i› power_mono ‹i ≥ 1›
        by (smt (verit, ccfv_SIG) of_nat_le_of_nat_power_cancel_iff)
      hence "N - a^3 = N - (int a)^3"
        by (simp add: of_nat_diff)
      thus ?thesis 
        using ‹odd (N - (int a)^3)› by presburger
    qed
    hence "odd (nat (N - a^3))"
      using nat_int by presburger
    moreover have "3*(nat k) ≥ 1"
      using ‹3 ≤ k› by auto
    ultimately have "∃b. odd b ∧ [(nat (N - a^3)) = b^3] (mod 2^(3*(nat k)))"
      using every_odd_nat_cong_cube by presburger
    then obtain b'::nat where "odd b'" and "[(nat (N - a^3)) = b'^3] (mod 2^(3*(nat k)))"
      by auto
    define b::nat where "b = b' mod (8^(nat k))"
    have "[nat (N - a^3) ≠ 0] (mod 2^(3*(nat k)))"
      by (meson ‹1 ≤ 3 * nat k› ‹odd (nat (int (N - a ^ 3)))› cong_0_iff even_power 
          cong_dvd_modulus_nat dvd_refl order_less_le_trans zero_less_one)
    hence "b > 0"
      by (metis ‹1 ≤ 3 * nat k› ‹2 ^ 3 = 8› b_def ‹odd b'› dvd_power gcd_nat.trans 
          mod_greater_zero_iff_not_dvd order_less_le_trans power_mult zero_less_one)
    hence "b^3 > 0"
      by simp

    have "b < 8 powr k"
      by (simp add: b_def powr_int)
    hence "b^3 < (8 powr (k))^3"
      using power_strict_mono by fastforce
    hence "b^3 < 8 powr (3*k)"
      by (metis mult.commute of_int_mult of_int_numeral powr_ge_zero powr_numeral powr_powr)
    have "N - a^3 ≥ 8*8 powr (3*k)"
      using ‹(i - 1)^3 < i^3› ‹8 * 8 powr (3*k) ≤ N - i^3›
      by (smt (verit, best)  a_def int_ops(6) of_int_less_iff of_int_of_nat_eq of_nat_power)
    
    have exp19:"7*8 powr (3*k) = 8*8 powr (3*k) - 8 powr (3*k)"
      by simp
    have exp20:"… < (N - a^3) - b^3"
      using ‹N - a^3 ≥ 8*8 powr (3*k)› ‹b^3 < 8 powr (3*k)› by linarith
    hence "… < N - a^3"
      using ‹0 < b ^ 3› ‹(b ^ 3) < 8 powr (3 * k)› by linarith
    have exp21: "… < 11*8 powr (3*k)"
      using a_def ‹(i-1)^3 < i^3› exp16 exp01 exp17 exp18 exp02
      by (smt(verit) of_int_less_iff of_int_of_nat_eq of_nat_diff of_nat_le_of_nat_power_cancel_iff 
          of_nat_power)

    have "[N - a^3 - b^3 = 0] (mod (8^(nat k)))"
      by (smt (verit, del_insts) ‹2^3 = 8› ‹[nat (N - a^3) = b'^3] (mod 2^(3*nat k))› b_def
          unique_euclidean_semiring_class.cong_def cong_diff_iff_cong_0_nat diff_is_0_eq' nat_int 
          nle_le power_mod power_mult)
    define q::real where "q = (N - a^3 - b^3)/((8 powr k))"
    have "(N - a^3 - b^3)/(8 powr k) ≥ 1"
      by (smt (verit, ccfv_SIG) ‹1 ≤ 3 * nat k› exp20 le_divide_eq_1_pos less_numeral_extra(1) 
          nat_0_less_mult_iff of_int_less_iff order_less_le_trans powr_gt_zero powr_less_cancel_iff 
          zero_less_nat_eq)
    hence "q ≠ 0"
      using q_def by auto
    moreover have "(N - a^3 - b^3) ≠ 0"
      using exp20 by auto

    moreover have "(8 powr k) ≠ 0"
      using power_not_zero by auto
    ultimately have exp22: "q*((8 powr k)) = (N - a^3 - b^3)"
      using q_def by auto

    have exp23: "(8::nat)^(nat k) = 8 powr k"
      using ‹1 ≤ 3 * nat k› powr_int by auto
    hence "q = ⌊q⌋"
      using ‹[N - a^3 - b^3 = 0] (mod (8::nat)^(nat k))› q_def
      by (metis cong_0_iff floor_of_nat of_int_of_nat_eq real_of_nat_div)
    have "7*8 powr (3*k) < q*8 powr (k)"
      using exp19 exp20 exp22
      by presburger
    hence "q > 7*8 powr (3*k)/(8 powr k)"
      by (simp add: divide_less_eq)
    hence "… > 7*(8 powr (3*k - k))"
      using powr_diff by (metis of_int_diff times_divide_eq_right)
    hence "… > 7 * 8 powr (2*k)"
      by simp

    have "q*8 powr k < 11*8 powr (3*k)"
      using exp22 exp21 by linarith
    hence "q < 11*8 powr (3*k) / 8 powr k"
      by (simp add: pos_less_divide_eq)
    hence "q < 11*8 powr (3*k - k)"
      using powr_diff by (metis of_int_diff times_divide_eq_right) 
    hence "q < 11*8 powr (2*k)"
      by simp
    have exp24:"(8::nat)^(2*(nat k)) = 8 powr (2*k)"
      using ‹3 ≤ k› powr_int
      by (metis exp23 of_int_mult of_int_numeral of_nat_power power_even_eq powr_ge_zero 
          powr_numeral powr_powr powr_powr_swap) 
    hence "6*(8::nat)^(2*(nat k)) = 6*8 powr (2*k)"
      by simp
    have "6*(8 powr (2* k)) < 7*(8 powr (2*k))"
      by simp
    hence "⌊q⌋ - 6*(8 ^ (2*(nat k))) > 0"
      using ‹7 * 8 powr (2 * k) < q› ‹q = ⌊q⌋›
    proof -
      have  "6 * 8 powr (2*k) < 7 * 8 powr (2*k)"
        by simp
      have "7 * 8 powr (2 * k) < q"
        using ‹7 * 8 powr (2 * k) < q› by fastforce
      hence "6 * 8 powr (2 * k) < q"
        using ‹6 * 8 powr (2*k) < 7 * 8 powr (2*k)› by linarith
      thus  ?thesis
        by (metis (no_types) ‹q = ⌊q⌋› exp24 of_int_0_less_iff of_int_diff of_int_eq_numeral_power_cancel_iff
            of_int_mult of_int_numeral real_of_nat_eq_numeral_power_cancel_iff diff_gt_0_iff_gt)
    qed
    then obtain r::nat where "r = ⌊q⌋ - 6*(8 ^ (2*(nat k)))"
      by (metis zero_less_imp_eq_int)
    hence "r = q - 6*(8 powr (2*k))"
      by (metis ‹q = ⌊q⌋› exp24 of_int_diff of_int_mult of_int_numeral of_int_of_nat_eq 
          real_of_nat_eq_numeral_power_cancel_iff)
    have "(22::nat)^3 < 8^6"
      using power_mono by simp
    have "… ≤ 8^(nat (k*2))"
      by (smt (verit, best) ‹k ≥ 3›  nat_mono nat_numeral numeral_Bit0 numeral_Bit1 
          numeral_eq_one_iff numeral_less_iff power_increasing_iff)
    have "… = 8 powr (2*k)"
      by (smt (verit, best) ‹3 ≤ k› numeral_Bit0 numeral_One of_nat_1 of_nat_numeral of_nat_power powr_int)
    have "… < r"
      using ‹r = q - 6*(8 powr (2*k))› ‹q > 7*8 powr (2*k)› by auto
    have "… < 5*(8 powr (2*k))"
      using ‹q < 11*8 powr (2*k)› ‹r = q - 6*(8 powr (2*k))› by auto

    have "r > 22^3"
      using ‹(22::nat)^3 < 8^6› ‹8^6 ≤ 8^(nat (k*2))› ‹real (8^(nat (k*2))) = 8 powr (2*k)› ‹8 powr (2*k) < r› 
      by linarith
    hence " r ≥ 10648"
      by simp
    then obtain d m where "d ≥ 0" and "d ≤ 22" and "r = d^3 + 6*m" and "is_sumpow 2 3 m"
      using values_geq_22_cubed_can_be_normalised by auto

    define A::nat where "A = 8^(nat k)"
    have "m ≤ r/6"
      using ‹r = d^3 + 6*m› ‹d ≥ 0› by linarith 
    have "… < (5*8^(nat (2*k)))/6"
      by (smt (verit, ccfv_SIG) ‹3 ≤ k› ‹r < 5*(8 powr (2*k))› divide_strict_right_mono powr_int)
    have "… < 8^(nat (2*k))"
      by simp
    have "… = A2"
      by (metis A_def ‹real (8^(2*nat k)) = 8 powr (2*k)› ‹real (8 ^ nat(k*2)) = 8 powr (2*k)› 
          mult.commute power_mult real_of_nat_eq_numeral_power_cancel_iff)

    define c::nat where "c = d*2^(nat k)"
    have "N = a^3 + b^3 + (8^(nat k))*q"
      by (smt (verit, del_insts) ‹1 ≤ 3 * nat k› ‹N - a^3 - b^3 ≠ 0› exp22
          diff_is_0_eq' mult.commute nat_0_iff nat_0_less_mult_iff nat_int nle_le of_nat_add 
          of_nat_diff powr_int zero_less_nat_eq zero_less_one)
    moreover have "… = a^3 + b^3 + 8^(nat k)*(6*8 powr (2*k) + r)"
      using ‹r = q - 6*(8 powr (2*k))› by simp
    moreover have "… = a^3 + b^3 + 8^(nat k)*(6*8^(nat(2*k)) + r)"
    proof -
      have "8 powr (int 2 * k) = (8 ^ nat (k * int 2))"
        using ‹real (8 ^ nat (k * 2)) = 8 powr (2 * k)› by auto
      thus ?thesis
        by simp
    qed
    moreover have "… = a^3 + b^3 + (8^(nat k)*d^3) + 8^(nat k)*(6*8^(nat(2*k)) + 6*m)"
      by (simp add: ‹r = d^3 + 6*m› add_mult_distrib2)
    moreover have "… = a^3 + b^3 + ((2^(nat k))*d)^3 + 8^(nat k)*(6*8^(nat(2*k)) + 6*m)"
      by (smt (verit) ‹2 ^ 3 = 8› mult.commute power_mult power_mult_distrib)
    moreover have "… = a^3 + b^3 + c^3 + 6*A*(A2 + m)"
      by (smt (z3) A_def c_def exp24 ‹real (8 ^ nat(k*2)) = 8 powr (2*k)› add_mult_distrib2
          mult.assoc mult.commute of_nat_eq_iff power_even_eq)
    ultimately have "N = a^3 + b^3 + c^3 + 6*A*(A2 + m)"
      using of_nat_eq_iff by metis
    have "is_sumpow 3 6 (6*A*(A2 + m))"
    proof -
      have "m ≤ A2"
        using ‹8 ^ nat (2*k) = real (A2)› ‹m ≤ r/6› ‹real r / 6 < 5*8^nat (2 * k)/6› by linarith
      thus ?thesis
        using sum_of_6_cubes ‹is_sumpow 2 3 m› by simp
    qed
    then obtain l::"nat list" where "length l = 6" and "6*A*(A2 + m) = sumpow 3 l"
      using is_sumpow_def by blast
    have "a^3 + b^3 + c^3 + sumpow 3 l = sumpow 3 (a#b#c#l)"
      by (simp add: fold_plus_sum_list_rev)
    hence "… = N"
      using ‹6 * A * (A2 + m) = sumpow 3 l› ‹N = a^3 + b^3 + c^3 + 6*A*(A2 + m)› by presburger
    have "length (a#b#c#l) = 9"
      using ‹length l = 6› by simp
    thus ?thesis
      using ‹sumpow 3 (a#b#c#l) = N› is_sumpow_def by blast
  next
    case leq8pow10
    consider (leq40000) "N ≤ 40000" | (ge40000) "N > 40000"
      by arith
    thus ?thesis
    proof (cases)
      case ge40000
      define a::int where "a = ⌊(N - 10000) powr (1/3)⌋"

      text‹The following inequalities differ from those in \cite{nathanson1996}, which erroneously result in
the false statement a > 31, we have corrected these mistakes.  This does not affect the rest of the proof.›
      have "(N - 10000) > 30000"
        using ge40000 by linarith
      hence "(N - 10000) powr (1/3) ≥ (30000) powr (1/3)"
        using powr_mono2 by simp
      hence "a ≥ ⌊(30000::int) powr (1/3)⌋"
        using a_def floor_mono by simp
      have "(30000::nat) > (31)^3"
        by simp
      hence "(30000::nat) powr (1/3) > (31^3) powr (1/3)"
        by (metis numeral_less_iff numeral_pow of_nat_numeral powr_less_mono2 zero_le_numeral 
            zero_less_divide_1_iff zero_less_numeral)
      hence "… > 31"
        by auto
      hence "⌊(30000::nat) powr (1/3)⌋ ≥ 31"
        by linarith
      hence "a ≥ 31"
        using ‹a ≥ ⌊(30000::int) powr (1/3)⌋›
        by (metis nle_le of_int_numeral of_nat_numeral order_trans)
      hence "nat a = a"
        by simp
      have "∀(x::int) > 4. x*(x - 3) > 1"
        by (simp add: less_1_mult)
      hence "∀(x::int) > 4. x^2 - 3*x -1 > 0"
        by (simp add: mult.commute power2_eq_square right_diff_distrib')
      define d::int where "d = (a+1)^3 - a^3"
      moreover have "… =3*a^2 + 3*a + 1"
        by Groebner_Basis.algebra
      moreover hence "a^2 - 3*a -1 > 0"
        using ‹a ≥ 31› ‹∀(x::int) > 4. x^2 - 3*x -1 > 0› by simp
      moreover hence "4*a^2 > 3*a^2 + 3*a + 1"
        by simp
      ultimately have "4*a^2 > d"
        by presburger
      have "a < (N powr (1/3))"
        by (smt (verit, best) a_def diff_less floor_eq_iff ge40000 of_nat_0_le_iff of_nat_less_iff 
            powr_less_mono2 zero_less_divide_1_iff zero_less_numeral)
      hence "a powr 2 < (N powr (1/3))powr 2"
        using ‹31 ≤ a› order_less_le by fastforce
      hence "4*a^2 < 4*(N powr (2/3))"
        using a_def powr_powr by fastforce
      have "N - (N - 10000) ≤ 10000"
        by simp
      hence "N - ((N - 10000) powr (1/3)) powr 3 ≤ 10000"
        using powr_powr powr_one_gt_zero_iff ge40000 by force
      hence "N - ⌊(N-10000) powr (1/3) + 1⌋ powr 3 < 10000"
        by (smt (verit, best) powr_ge_zero powr_less_mono2 real_of_int_floor_add_one_gt)
      hence "N - (a+1) powr 3 < 10000"
        using a_def one_add_floor by simp
      have "(a+1) powr 3 = (a+1)^3"
        using a_def by auto
      hence "N - (a+1)^3 < 10000"
        using ‹N - (a+1) powr 3 < 10000› by linarith
      have "… ≤ N - (N - 10000)"
        using ge40000 by auto 
      moreover have "… = N - ((N - 10000) powr (1/3)) powr 3"
        using powr_powr by simp
      moreover have "… ≤ N - ⌊((N - 10000) powr (1/3))⌋ powr 3"
        by simp
      moreover have "… = N - a^3"
        using a_def by simp
      ultimately have "10000 ≤ N - a^3"
        by linarith
      have "… = N - (a+1)^3 + d"
        using d_def by simp
      moreover have "… < 10000 + 4*(N powr (2/3))"
        using ‹N - (a+1)^3 < 10000› ‹4*a^2 > d› ‹4*a^2 < 4*(N powr (2/3))› le_less_trans by linarith
      ultimately have "N - a^3 < 10000 + 4*(N powr (2/3))"
        by presburger
      hence "int (N - (nat a)^3) < 10000 + 4*(N powr (2/3))"
        using ‹nat a = a›
        by (smt (verit, best) ‹10000 ≤ int N - a^3› int_ops(6) of_int_of_nat_eq of_nat_power)

      consider (N_min_a_cube_leq40000) "N-(nat a)^3 ≤ 40000" | (N_min_a_cube_ge40000) "N - (nat a)^3 > 40000"
        by arith
      thus ?thesis
      proof (cases)
        case N_min_a_cube_leq40000
        have "N - (nat a)^3 ≥ 10000"
          using ‹10000 ≤ N - a^3› int_nat_eq ‹a ≥ 31›
          by (smt (verit) int_ops(6) numeral_Bit0 numeral_Bit1 numeral_One of_nat_1 of_nat_le_iff 
              of_nat_numeral of_nat_power) 
        hence "is_sumpow 3 6 (N-(nat a)^3)"
          using leq_40000_is_sum_of_9_cubes N_min_a_cube_leq40000 by force

        then obtain l::"nat list" where "length l = 6" and "N - (nat a)^3 = sumpow 3 l"
          using is_sumpow_def by blast
        have "0 + 0 + (nat a)^3 + sumpow 3 l = sumpow 3 ((nat a)#0#0#l)"
          by (simp add: fold_plus_sum_list_rev)
        hence "… = N"
          using ‹N - (nat a)^3 = sumpow 3 l› ‹10000 ≤ N - nat a ^ 3› by presburger 
        have "length ((nat a)#0#0#l) = 9"
          using ‹length l = 6› by simp
        thus ?thesis
          using ‹sumpow 3 ((nat a)#0#0#l) = N› is_sumpow_def by blast
      next
        case N_min_a_cube_ge40000
        define N'::nat where "N' = N - (nat a)^3"
        hence "N' > 40000"
          using N_min_a_cube_ge40000 by simp

        define b::int where "b = ⌊(N' - 10000) powr (1/3)⌋"
        text‹The same mistake as above crops up, and we have corrected it in the same way.›
        have "… ≥ 31"
          by (smt (verit) ‹31 ≤ ⌊real 30000 powr (1 / 3)⌋› ‹40000 < N'› floor_mono floor_of_nat 
              int_ops(2) int_ops(6) less_imp_of_nat_less numeral_Bit0 numeral_Bit1 numeral_One 
              of_nat_numeral powr_mono2 zero_le_divide_1_iff)
        hence "b ≥ 31"
          using b_def by simp
        hence "nat b = b"
          by simp

        define d'::int where "d' = (b+1)^3 - b^3"
        moreover have "… =3*b^2 + 3*b + 1"
          by Groebner_Basis.algebra 
        moreover have "b^2 - 3*b -1 > 0"
          using ‹b ≥ 31› ‹∀(x::int) > 4. x^2 - 3*x -1 > 0› by simp
        moreover hence "4*b^2 > 3*b^2 + 3*b + 1"
          by simp
        ultimately have "4*b^2 > d'"        
          by presburger

        have "b < (N' powr (1/3))"
          using b_def
          by (smt (verit, best) N_min_a_cube_ge40000 N'_def diff_less floor_eq_iff of_nat_0_le_iff
              of_nat_less_iff powr_less_mono2 zero_less_divide_1_iff zero_less_numeral)
        hence "b powr 2 < (N' powr (1/3))powr 2"
          using ‹31 ≤ b› order_less_le by fastforce
        hence "4*b^2 < 4*(N' powr (2/3))"
          using b_def powr_powr by fastforce
        have "N' - (N' - 10000) ≤ 10000"
          by simp
        hence "N' - ((N' - 10000) powr (1/3)) powr 3 ≤ 10000"
          using powr_powr powr_one_gt_zero_iff ‹N' > 40000› by force
        hence "N' - ⌊(N'-10000) powr (1/3) + 1⌋ powr 3 < 10000"
          by (smt (verit, best) powr_ge_zero powr_less_mono2 real_of_int_floor_add_one_gt)
        hence "N' - (b+1) powr 3 < 10000"
          using b_def one_add_floor by simp
        have "(b+1) powr 3 = (b+1)^3"
          using b_def by auto
        hence "N' - (b+1)^3 < 10000"
          using  ‹N' - (b+1) powr 3 < 10000› by linarith
        have "… ≤ N' - (N' - 10000)"
          using ‹N' > 40000› by auto
        moreover have "… = N' - ( (N' - 10000) powr (1/3)) powr 3"
          using powr_powr by simp
        moreover have "… ≤ N' - ⌊((N' - 10000) powr (1/3))⌋ powr 3"
          by simp
        moreover have "… = N' - b^3"
          using b_def by simp
        ultimately have "10000 ≤ N' - b^3"
          by linarith
        have "… = N' - (b+1)^3 + d'"
          using d'_def by simp
        moreover have "… < 10000 + 4*(N' powr (2/3))"
          using ‹N' - (b+1)^3 < 10000› ‹4*b^2 > d'› ‹4*b^2 < 4*(N' powr (2/3))› le_less_trans 
          by linarith
        ultimately have "N' - b^3 < 10000 + 4*(N' powr (2/3))"
          by presburger
        hence "int (N' - (nat b)^3) < 10000 + 4*(N' powr (2/3))"
          using ‹nat b = b›
          by (smt (verit, best) ‹10000 ≤ int N' - b ^ 3› int_ops(6) of_int_of_nat_eq of_nat_power)

        consider (N'_min_b_cube_leq40000) "N'-(nat b)^3 ≤ 40000" | (N'_min_b_cube_ge40000) "N' - (nat b)^3 > 40000"
          by arith
        thus ?thesis
        proof (cases)
          case N'_min_b_cube_leq40000
          have "N' - (nat b)^3 ≥ 10000"
            using ‹10000 ≤ N' - b^3› int_nat_eq ‹b ≥ 31›
            by (smt (verit) int_ops(6) numeral_Bit0 numeral_Bit1 numeral_One of_nat_1 of_nat_le_iff 
                of_nat_numeral of_nat_power) 
          hence "N - (nat a)^3 - (nat b)^3 ≥ 10000"
            using N'_def by simp
          hence "is_sumpow 3 6 (N - (nat a)^3 - (nat b)^3)"
            using leq_40000_is_sum_of_9_cubes N'_min_b_cube_leq40000 N'_def by force
          then obtain l::"nat list" where "length l = 6" and "N - (nat a)^3 - (nat b)^3 = sumpow 3 l"
            using is_sumpow_def by blast
          have "0 + (nat b)^3 + (nat a)^3 + sumpow 3 l = sumpow 3 ((nat a)#(nat b)#0#l)"
            by (simp add: fold_plus_sum_list_rev)
          hence "… = N"
            using ‹N - (nat a)^3 - (nat b)^3 = sumpow 3 l› ‹10000 ≤ N - nat a^3 - nat b^3› 
            by presburger 
          have "length ((nat a)#(nat b)#0#l) = 9"
            using ‹length l = 6› by simp
          thus ?thesis 
            using ‹sumpow 3 (nat a # nat b # 0 # l) = N› is_sumpow_def by blast
        next
          case N'_min_b_cube_ge40000

          define N''::nat where "N'' = N - (nat a)^3 - (nat b)^3"
          hence "N'' > 40000"
            using N'_min_b_cube_ge40000 N'_def by simp

          define c::int where "c = ⌊(N'' - 10000) powr (1/3)⌋"
          text‹We correct the same mistake as above.›
          have "… ≥ 31"
            by (smt (verit) ‹31 ≤ ⌊real 30000 powr (1 / 3)⌋› ‹40000 < N''› floor_mono floor_of_nat 
                int_ops(2) int_ops(6) less_imp_of_nat_less numeral_Bit0 numeral_Bit1 numeral_One 
                of_nat_numeral powr_mono2 zero_le_divide_1_iff)
          hence "c ≥ 31"
            using c_def by simp
          hence cnotneg: "nat c = c"
            using int_nat_eq by simp 

          define d''::int where "d'' = (c+1)^3 - c^3"
          moreover have "… =3*c^2 + 3*c + 1"
            by Groebner_Basis.algebra 
          moreover have "c^2 - 3*c -1 > 0"
            using ‹c ≥ 31› ‹∀(x::int) > 4. x^2 - 3*x -1 > 0› by simp
          moreover hence "4*c^2 > 3*c^2 + 3*c + 1"
            by simp
          ultimately have "4*c^2 > d''"        
            by presburger
      
          have "c < (N'' powr (1/3))"
            using c_def N'_min_b_cube_ge40000 N'_def N''_def
            by (smt (verit, ccfv_SIG) diff_less floor_eq_iff of_nat_0_le_iff of_nat_less_iff 
                powr_less_mono2 zero_less_divide_1_iff zero_less_numeral)
          hence "c powr 2 < (N'' powr (1/3))powr 2"
            using ‹31 ≤ c› order_less_le by fastforce
          hence "4*c^2 < 4*(N'' powr (2/3))"
            using c_def powr_powr by fastforce

          have "N'' - (N'' - 10000) ≤ 10000"
            by simp
          hence "N'' - ((N'' - 10000) powr (1/3)) powr 3 ≤ 10000"
            using powr_powr powr_one_gt_zero_iff ‹N' > 40000› by force
          hence "N'' - ⌊(N''-10000) powr (1/3) + 1⌋ powr 3 < 10000"
            by (smt (verit, best) powr_ge_zero powr_less_mono2 real_of_int_floor_add_one_gt)
          hence "N'' - (c+1) powr 3 < 10000"
            using c_def one_add_floor by simp
          have "(c+1) powr 3 = (c+1)^3"
            using c_def by auto
          hence "N'' - (c+1)^3 < 10000"
            using  ‹N'' - (c+1) powr 3 < 10000› by linarith
          have "… ≤ N'' - (N'' - 10000)"
            using ‹N'' > 40000› by auto
          moreover have "… = N'' - ( (N'' - 10000) powr (1/3)) powr 3"
            using powr_powr by simp
          moreover have "… ≤ N'' - ⌊((N'' - 10000) powr (1/3))⌋ powr 3"
            by simp
          moreover have "… = N'' - c^3"
            using c_def by simp
          ultimately have "10000 ≤ N'' - c^3"
            by linarith
          moreover have "… = N'' - (c+1)^3 + d''"
            using d''_def by simp
          moreover have "… < 10000 + 4*(N'' powr (2/3))"
            using ‹N'' - (c+1)^3 < 10000› ‹4*c^2 > d''› ‹4*c^2 < 4*(N'' powr (2/3))› le_less_trans 
            by linarith
          ultimately have "N'' - c^3 < 10000 + 4*(N'' powr (2/3))"
            by presburger
          hence "int (N'' - (nat c)^3) < 10000 + 4*(N'' powr (2/3))"
            using cnotneg
            by (smt (verit, best) ‹10000 ≤ int N'' - c ^ 3› int_ops(6) of_int_of_nat_eq of_nat_power) 
      
          have "N - (nat a)^3 - (nat b)^3 - (c+1)^3 < 10000"
            using ‹N'' - (c+1)^3 < 10000› N''_def by simp
          have "… ≤ N'' - (nat c)^3"
            by (smt (verit, del_insts) ‹10000 ≤ int N'' - c ^ 3› cnotneg int_ops(6) of_nat_power)
          have "… < 10000 + 4*((N - (nat a)^3 - (nat b)^3) powr (2/3))"
            using N''_def ‹int (N'' - (nat c)^3) < 10000 + 4*(N'' powr (2/3))› by simp
          moreover have "… < 10000 + 4*((10000 + 4*((N - (nat a)^3) powr (2/3))) powr (2/3))" 
            using ‹int (N' - (nat b)^3) < 10000 + 4*(N' powr (2/3))› powr_less_mono2 N'_def by simp
          moreover have "… < 10000 + 4*((10000 + 4*((10000 + 4*(N powr (2/3))) powr (2/3))) powr (2/3))"
            using ‹int (N - (nat a)^3) < 10000 + 4*(N powr (2/3))› powr_less_mono2 by simp
          moreover have "… ≤ 10000 + 4*((10000 + 4*((10000 + 4*((int 8^10) powr (2/3))) powr (2/3))) powr (2/3))"
            using leq8pow10 powr_less_mono2 nle_le
            by (smt (verit, best) numeral_power_le_of_nat_cancel_iff of_int_of_nat_eq of_nat_0_le_iff 
                of_nat_power powr_ge_zero zero_less_divide_iff)
          moreover have "… = 10000 + 4*((10000 + 4*((4204304) powr (2/3))) powr (2/3))"
            using powr_powr by simp
          (*4251528 is the next lowest cube number after 4204304*)
          moreover have "… ≤ 10000 + 4*((10000 + 4*((4251528) powr (2/3))) powr (2/3))" 
            by (smt (verit, best) powr_ge_zero powr_less_mono2 zero_less_divide_iff)
          moreover have "… = 10000 + 4*((114976) powr (2/3))"
            by auto
          (*117649 is the next lowest cube number after 114976*)
          moreover have "… ≤ 10000 + 4*((117649) powr (2/3))" 
            by (smt (verit, best) powr_ge_zero powr_less_mono2 zero_less_divide_iff)
          moreover have "… ≤ 20000"
            by auto
          ultimately have "N - (nat a)^3 - (nat b)^3 - (nat c)^3 ≤ 20000"
            by (smt (verit, del_insts) N''_def numeral_Bit0 numeral_Bit1 numerals(1) of_int_of_nat_eq 
                of_nat_1 of_nat_le_iff of_nat_numeral)

          hence "is_sumpow 3 6 (N - (nat a)^3 - (nat b)^3 - (nat c)^3)"
            using leq_40000_is_sum_of_9_cubes ‹10000 ≤ N'' - c^3› N''_def cnotneg 
            by (smt (verit) ‹10000 ≤ int (N'' - nat c ^ 3)› int_ops(2) numeral_Bit0 numeral_Bit1 
                numeral_One of_nat_le_iff of_nat_numeral order_less_le) 
          then obtain l::"nat list" where "length l = 6" and 
            "N - (nat a)^3 - (nat b)^3 - (nat c)^3 = sumpow 3 l"
            using is_sumpow_def by blast
          moreover have "(nat c)^3 + (nat b)^3 + (nat a)^3 + sumpow 3 l = sumpow 3 ((nat a)#(nat b)#(nat c)#l)"
            by (simp add: fold_plus_sum_list_rev)
          ultimately have  "… = N"
            using ‹10000 ≤ int N'' - c^3› N''_def cnotneg 
            by (smt (verit, del_insts) diff_diff_left of_nat_power of_nat_le_iff le_add_diff_inverse 
                add.commute int_ops(6))
          moreover have "length ((nat a)#(nat b)#(nat c)#l) = 9"
            using ‹length l = 6› by simp
          ultimately show ?thesis 
            using is_sumpow_def by blast
        qed
      qed
    next
      case leq40000
      thus ?thesis 
        using leq_40000_is_sum_of_9_cubes by blast
    qed
  qed
  thus ?thesis
    by simp
qed
end