2607.19283v1 / ENOTV/EvenEnergy.lean

all files

import ENOTV.Junction
import ENOTV.SourceComparison

/-!
# Local structure of the alternating Euler construction

This file proves the pointwise block and junction descriptions needed for
the global even-order energy estimate.
-/

noncomputable section

open scoped BigOperators

namespace ENOTV

theorem translated_pairJumps_zero_before (d L p : ℕ) (i : ℤ)
    (hi : i < (4 * L * p : ℕ)) :
    translate (4 * L * p) (pairJumps d L) i = 0 := by
  by_contra h
  have hs : i ∈ (translate (4 * L * p) (pairJumps d L)).support :=
    Finsupp.mem_support_iff.mpr h
  have hs' : i - (4 * L * p : ℕ) ∈ (pairJumps d L).support := by
    simpa [Finsupp.mem_support_iff] using hs
  have hb := pairJumps_support d L hs'
  simp only [Finset.mem_Ico] at hb
  omega

theorem translated_pairJumps_zero_after (d L p : ℕ) (i : ℤ)
    (hi : (4 * L * (p + 1) : ℕ) ≤ i) :
    translate (4 * L * p) (pairJumps d L) i = 0 := by
  by_contra h
  have hs : i ∈ (translate (4 * L * p) (pairJumps d L)).support :=
    Finsupp.mem_support_iff.mpr h
  have hs' : i - (4 * L * p : ℕ) ∈ (pairJumps d L).support := by
    simpa [Finsupp.mem_support_iff] using hs
  have hb := pairJumps_support d L hs'
  simp only [Finset.mem_Ico] at hb
  push_cast at hi hb
  ring_nf at hi hb
  omega

theorem multiJumps_eq_pair_on_interval {d L P p : ℕ} (hL : 0 < L)
    (hp : p < P) (i : ℤ)
    (hi0 : (4 * L * p : ℕ) ≤ i)
    (hi1 : i < (4 * L * (p + 1) : ℕ)) :
    multiJumps d L P i =
      translate (4 * L * p) (pairJumps d L) i := by
  simp only [multiJumps, Finsupp.finsetSum_apply]
  apply Finset.sum_eq_single p
  · intro q hq hqp
    simp only [Finset.mem_range] at hq
    by_cases hqp' : q < p
    · apply translated_pairJumps_zero_after
      have hnat : 4 * L * (q + 1) ≤ 4 * L * p := by
        apply Nat.mul_le_mul_left
        omega
      have hcast : (4 * L * (q + 1) : ℕ) ≤
          (4 * L * p : ℕ) := hnat
      have hcast' : ((4 * L * (q + 1) : ℕ) : ℤ) ≤
          ((4 * L * p : ℕ) : ℤ) := by exact_mod_cast hcast
      exact hcast'.trans hi0
    · have hpq : p < q := by omega
      apply translated_pairJumps_zero_before
      have hnat : 4 * L * (p + 1) ≤ 4 * L * q := by
        apply Nat.mul_le_mul_left
        omega
      have hsep : i < (4 * L * q : ℕ) := by
        exact hi1.trans_le (by exact_mod_cast hnat)
      exact hsep
  · intro hpnot
    exact (hpnot (by simpa using hp)).elim

theorem pairJumps_first_block (d L t : ℕ) (ht : t < 2 * L) :
    pairJumps d L t = realBlock d L t := by
  simp [pairJumps, pairJumpNat_apply, ht]

theorem pairJumps_second_block (d L t : ℕ) (hL : 0 < L)
    (ht : t < 2 * L) :
    pairJumps d L (2 * L + t) = -realBlock d L t := by
  have h2 : ¬2 * L + t < 2 * L := by omega
  have h4 : 2 * L + t < 4 * L := by omega
  change liftNat (pairJumpNat d L) (2 * (L : ℤ) + (t : ℤ)) =
    -realBlock d L t
  rw [show 2 * (L : ℤ) + (t : ℤ) =
    ((2 * L + t : ℕ) : ℤ) by push_cast; ring,
    liftNat_apply_nat, pairJumpNat_apply]
  simp [h2, h4]

/-- Exact value on any of the `2P` alternating blocks. -/
theorem multiJumps_block_value {d L P n t : ℕ} (hL : 0 < L)
    (hn : n < 2 * P) (ht : t < 2 * L) :
    multiJumps d L P (2 * L * n + t : ℕ) =
      (-1 : ℝ) ^ n * realBlock d L t := by
  rcases Nat.even_or_odd n with hnEven | hnOdd
  · obtain ⟨p, hp⟩ := hnEven
    rw [hp]
    have hpP : p < P := by omega
    rw [multiJumps_eq_pair_on_interval hL hpP]
    · rw [translate_apply]
      have harg :
          ((2 * L * (p + p) + t : ℕ) : ℤ) -
              4 * (L : ℤ) * (p : ℤ) = (t : ℤ) := by
        push_cast
        ring
      rw [harg]
      have hpone : (-1 : ℝ) ^ (p + p) = 1 := by
        rw [pow_add, ← mul_pow]
        norm_num
      rw [hpone, one_mul]
      exact pairJumps_first_block d L t ht
    · have ht0 : (0 : ℤ) ≤ t := by positivity
      push_cast
      ring_nf
      nlinarith
    · have ht' : (t : ℤ) < 2 * L := by exact_mod_cast ht
      push_cast
      ring_nf
      nlinarith
  · obtain ⟨p, hp⟩ := hnOdd
    rw [hp]
    have hpP : p < P := by omega
    rw [multiJumps_eq_pair_on_interval hL hpP]
    · have harg :
          ((2 * L * (2 * p + 1) + t : ℕ) : ℤ) -
              4 * (L : ℤ) * (p : ℤ) =
            2 * (L : ℤ) + (t : ℤ) := by
        push_cast
        ring
      rw [translate_apply, harg]
      simpa [pow_add, pow_succ] using pairJumps_second_block d L t hL ht
    · have ht0 : (0 : ℤ) ≤ t := by positivity
      push_cast
      ring_nf
      nlinarith
    · have ht' : (t : ℤ) < 2 * L := by exact_mod_cast ht
      push_cast
      ring_nf
      nlinarith

theorem multiJumps_zero_of_neg (d L P : ℕ) {i : ℤ} (hi : i < 0) :
    multiJumps d L P i = 0 := by
  simp only [multiJumps, Finsupp.finsetSum_apply]
  apply Finset.sum_eq_zero
  intro p hp
  apply translated_pairJumps_zero_before
  have hnonneg : (0 : ℤ) ≤ (4 * L * p : ℕ) := by positivity
  omega

theorem multiJumps_zero_of_ge (d L P : ℕ) {i : ℤ}
    (hi : (4 * L * P : ℕ) ≤ i) :
    multiJumps d L P i = 0 := by
  simp only [multiJumps, Finsupp.finsetSum_apply]
  apply Finset.sum_eq_zero
  intro p hp
  simp only [Finset.mem_range] at hp
  apply translated_pairJumps_zero_after
  have hnat : 4 * L * (p + 1) ≤ 4 * L * P := by
    apply Nat.mul_le_mul_left
    omega
  have hcast : ((4 * L * (p + 1) : ℕ) : ℤ) ≤
      ((4 * L * P : ℕ) : ℤ) := by exact_mod_cast hnat
  exact hcast.trans hi

/-! ## Junction coordinates -/

def realJunctionProfile (d L : ℕ) (h : ℤ) : ℝ :=
  (junctionProfile d L (h : ℚ) : ℚ)

theorem diffIter_congr_of_block (n : ℕ) (f g : Seq) (j : ℤ)
    (h : ∀ q : ℕ, q ≤ n → f (j + q) = g (j + q)) :
    diffIter n f j = diffIter n g j := by
  rw [diffIter_eq_fwdDiff, diffIter_eq_fwdDiff,
    fwdDiff_iter_eq_sum_shift, fwdDiff_iter_eq_sum_shift]
  apply Finset.sum_congr rfl
  intro q hq
  simp only [Finset.mem_range] at hq
  congr 1
  simpa using h q (by omega)

theorem diffIter_zero_seq (n : ℕ) (j : ℤ) :
    diffIter n (0 : Seq) j = 0 := by
  induction n generalizing j with
  | zero => rfl
  | succ n ih =>
      rw [diffIter, diff, ih, ih, sub_zero]

theorem diffIter_mul_comp_sub (n : ℕ) (c : ℝ) (g : Seq) (B s : ℤ) :
    diffIter n (fun x => c * g (x - B)) (B + s) =
      c * diffIter n g s := by
  rw [diffIter_eq_fwdDiff, diffIter_eq_fwdDiff,
    fwdDiff_iter_eq_sum_shift, fwdDiff_iter_eq_sum_shift,
    Finset.mul_sum]
  apply Finset.sum_congr rfl
  intro q hq
  simp only [zsmul_eq_mul, nsmul_eq_mul, mul_one]
  ring

theorem realJunction_high_difference (d L : ℕ) (s : ℤ) :
    diffIter (d + 1) (realJunctionProfile d L) s =
      ((fwdDiff (1 : ℚ))^[d + 1] (junctionProfile d L) (s : ℚ) : ℚ) := by
  rw [diffIter_eq_fwdDiff, fwdDiff_iter_eq_sum_shift,
    fwdDiff_iter_eq_sum_shift]
  unfold realJunctionProfile
  push_cast
  apply Finset.sum_congr rfl
  intro q hq
  simp only [zsmul_eq_mul, nsmul_eq_mul, mul_one]
  norm_cast

theorem multiJumps_junction_value {d L P n : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : 0 < L)
    (hn0 : 1 ≤ n) (hnP : n < 2 * P) (h : ℤ)
    (hlo : -(2 * (L : ℤ)) ≤ h) (hhi : h < 2 * (L : ℤ)) :
    multiJumps d L P ((2 * L * n : ℕ) + h) =
      (-1 : ℝ) ^ (n - 1) * realJunctionProfile d L h := by
  by_cases hhneg : h < 0
  · let t := Int.toNat (2 * (L : ℤ) + h)
    have ht0 : 0 ≤ 2 * (L : ℤ) + h := by omega
    have htcast : (t : ℤ) = 2 * (L : ℤ) + h :=
      Int.toNat_of_nonneg ht0
    have ht : t < 2 * L := by
      rw [Int.toNat_lt ht0]
      push_cast
      omega
    have hnprev : n - 1 < 2 * P := by omega
    have hindex :
        ((2 * L * n : ℕ) : ℤ) + h =
          ((2 * L * (n - 1) + t : ℕ) : ℤ) := by
      push_cast
      rw [Nat.cast_sub hn0]
      rw [htcast]
      ring
    rw [hindex, multiJumps_block_value hL hnprev ht]
    have hblock :
        realBlock d L t =
          (blockPoly d L (2 * (L : ℤ) + h : ℤ) : ℚ) := by
      unfold realBlock
      have hq : (t : ℚ) = (2 * (L : ℤ) + h : ℤ) := by
        exact_mod_cast htcast
      rw [hq]
    rw [hblock]
    simp [realJunctionProfile, junctionProfile, le_of_lt hhneg]
  · by_cases hh0 : h = 0
    · subst h
      have hn : n < 2 * P := hnP
      have hz := multiJumps_block_value (d := d) hL hn
        (show 0 < 2 * L by omega)
      rw [show ((2 * L * n : ℕ) : ℤ) + 0 =
          ((2 * L * n + 0 : ℕ) : ℤ) by simp, hz]
      simp [realBlock, blockPoly_zero d L hd heven,
        realJunctionProfile, junctionProfile,
        blockPoly_two_mul d L hd heven]
    · have hhpos : 0 < h := by omega
      let t := Int.toNat h
      have htcast : (t : ℤ) = h := Int.toNat_of_nonneg hhpos.le
      have ht : t < 2 * L := by
        rw [Int.toNat_lt hhpos.le]
        exact hhi
      have hindex :
          ((2 * L * n : ℕ) : ℤ) + h =
            ((2 * L * n + t : ℕ) : ℤ) := by
        push_cast
        rw [htcast]
      rw [hindex, multiJumps_block_value hL hnP ht]
      have hblock :
          realBlock d L t = (blockPoly d L (h : ℤ) : ℚ) := by
        unfold realBlock
        have hq : (t : ℚ) = (h : ℚ) := by
          exact_mod_cast htcast
        rw [hq]
      rw [hblock]
      have hnform : n = (n - 1) + 1 := by omega
      rw [hnform, pow_succ]
      simp [realJunctionProfile, junctionProfile, not_le_of_gt hhpos]

theorem multiJumps_internal_difference_le {d L P n : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hn0 : 1 ≤ n) (hnP : n < 2 * P) (s : ℤ)
    (hs0 : -(d + 1 : ℕ) ≤ s) (hs1 : s ≤ -1) :
    |cdiffIter (d + 1) (multiJumps d L P)
        ((2 * L * n : ℕ) + s)| ≤ (junctionUpper d : ℝ) := by
  have hLpos : 0 < L := by omega
  rw [cdiffIter_apply]
  let B : ℤ := (2 * L * n : ℕ)
  let c : ℝ := (-1 : ℝ) ^ (n - 1)
  let g : Seq := fun x => c * realJunctionProfile d L (x - B)
  have hcongr :
      diffIter (d + 1) (multiJumps d L P : Seq) (B + s) =
        diffIter (d + 1) g (B + s) := by
    apply diffIter_congr_of_block
    intro q hq
    have hq' : q < d + 2 := by omega
    have hlo :
        -(2 * (L : ℤ)) ≤ s + (q : ℤ) := by
      push_cast
      omega
    have hhi :
        s + (q : ℤ) < 2 * (L : ℤ) := by
      push_cast
      omega
    have hv := multiJumps_junction_value hd heven hLpos hn0 hnP
      (s + (q : ℤ)) hlo hhi
    change multiJumps d L P (B + s + q) = g (B + s + q)
    rw [show B + s + (q : ℤ) = B + (s + (q : ℤ)) by ring, hv]
    simp [g, c, B]
  rw [hcongr, diffIter_mul_comp_sub]
  have hr := junction_difference_le d L hLpos s
    (by omega) (by omega)
  rw [realJunction_high_difference]
  have hcast :
      |(((fwdDiff (1 : ℚ))^[d + 1] (junctionProfile d L)
          (s : ℚ) : ℚ) : ℝ)| ≤ (junctionUpper d : ℝ) := by
    exact_mod_cast hr
  simpa [c] using hcast

/-! ## Amplitudes near junctions and outer endpoints -/

def boundaryAmplitudeUpper (d L : ℕ) : ℝ :=
  (endpointUpper d (d + 1) : ℝ) * (L : ℝ) ^ (d - 1)

theorem abs_realJunctionProfile_le {d L : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (h : ℤ) (hh0 : -(d + 1 : ℕ) ≤ h) (hh1 : h ≤ d + 1) :
    |realJunctionProfile d L h| ≤ boundaryAmplitudeUpper d L := by
  have hLpos : 0 < L := by omega
  by_cases hh : h ≤ 0
  · let r := Int.toNat (-h)
    have hr0 : 0 ≤ -h := by omega
    have hrcast : (r : ℤ) = -h := Int.toNat_of_nonneg hr0
    have hrR : r ≤ d + 1 := by
      rw [← Int.ofNat_le]
      push_cast
      omega
    have hrinside : r ≤ 2 * L := by omega
    have hb := abs_blockPoly_right_endpoint_le d (d + 1) L r
      hd heven hLpos hrR hrinside
    have hrcastQ : (r : ℚ) = -(h : ℚ) := by
      exact_mod_cast hrcast
    have harg : 2 * (L : ℚ) + (h : ℚ) =
        2 * (L : ℚ) - (r : ℚ) := by
      rw [hrcastQ]
      ring
    unfold realJunctionProfile junctionProfile boundaryAmplitudeUpper
    rw [if_pos (by exact_mod_cast hh)]
    rw [harg]
    exact_mod_cast hb
  · have hhpos : 0 < h := by omega
    let r := Int.toNat h
    have hrcast : (r : ℤ) = h := Int.toNat_of_nonneg hhpos.le
    have hrR : r ≤ d + 1 := by
      rw [← Int.ofNat_le]
      push_cast
      omega
    have hrinside : r ≤ 2 * L := by omega
    have hb := abs_blockPoly_endpoint_le d (d + 1) L r
      hd heven hLpos hrR hrinside
    unfold realJunctionProfile junctionProfile boundaryAmplitudeUpper
    rw [if_neg (by exact_mod_cast (not_le_of_gt hhpos))]
    have harg : (r : ℚ) = (h : ℚ) := by exact_mod_cast hrcast
    have hb' :
        |((blockPoly d L (h : ℚ) : ℚ) : ℝ)| ≤
          (endpointUpper d (d + 1) : ℝ) * (L : ℝ) ^ (d - 1) := by
      rw [← harg]
      exact_mod_cast hb
    simpa using hb'

theorem internal_blockAmplitude_le {d L P n : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hn0 : 1 ≤ n) (hnP : n < 2 * P) (s : ℤ)
    (hs0 : -(d + 1 : ℕ) ≤ s) (hs1 : s ≤ -1) :
    blockAmplitude (d + 2) (multiJumps d L P : Seq)
      ((2 * L * n : ℕ) + s) ≤ boundaryAmplitudeUpper d L := by
  unfold blockAmplitude
  simp only [dif_pos (by omega : 0 < d + 2)]
  apply Finset.sup'_le
  intro q hq
  simp only [Finset.mem_range] at hq
  have hlo : -(2 * (L : ℤ)) ≤ s + (q : ℤ) := by
    push_cast
    omega
  have hhi : s + (q : ℤ) < 2 * (L : ℤ) := by
    push_cast
    omega
  have hv := multiJumps_junction_value hd heven (by omega : 0 < L)
    hn0 hnP (s + (q : ℤ)) hlo hhi
  rw [show ((2 * L * n : ℕ) : ℤ) + s + (q : ℤ) =
      ((2 * L * n : ℕ) : ℤ) + (s + (q : ℤ)) by ring, hv,
    abs_mul, abs_neg_one_pow, one_mul]
  apply abs_realJunctionProfile_le hd heven hL
  · push_cast
    omega
  · push_cast
    omega

theorem multiJumps_left_endpoint_value_le {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hP : 0 < P) (h : ℤ) (hh0 : -(d + 1 : ℕ) ≤ h)
    (hh1 : h ≤ d + 1) :
    |multiJumps d L P h| ≤ boundaryAmplitudeUpper d L := by
  by_cases hh : h < 0
  · rw [multiJumps_zero_of_neg d L P hh, abs_zero]
    unfold boundaryAmplitudeUpper
    exact mul_nonneg (by exact_mod_cast endpointUpper_nonneg d (d + 1))
      (by positivity)
  · let t := Int.toNat h
    have htcast : (t : ℤ) = h := Int.toNat_of_nonneg (by omega)
    have ht : t < 2 * L := by
      rw [Int.toNat_lt (by omega)]
      push_cast
      omega
    have hn : 0 < 2 * P := by omega
    have hv := multiJumps_block_value (d := d) (n := 0) (t := t)
      (by omega : 0 < L) hn ht
    have hv' : multiJumps d L P (t : ℤ) = realBlock d L t := by
      simpa using hv
    rw [← htcast, hv']
    exact abs_realBlock_endpoint_le d (d + 1) L t hd heven
      (by omega) (by
        rw [← Int.ofNat_le]
        push_cast
        omega) (by omega)

theorem multiJumps_right_endpoint_value_le {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hP : 0 < P) (h : ℤ) (hh0 : -(d + 1 : ℕ) ≤ h)
    (hh1 : h ≤ d + 1) :
    |multiJumps d L P ((4 * L * P : ℕ) + h)| ≤
      boundaryAmplitudeUpper d L := by
  by_cases hh : 0 ≤ h
  · rw [multiJumps_zero_of_ge d L P (by
      push_cast
      omega), abs_zero]
    unfold boundaryAmplitudeUpper
    exact mul_nonneg (by exact_mod_cast endpointUpper_nonneg d (d + 1))
      (by positivity)
  · let r := Int.toNat (-h)
    let t := 2 * L - r
    have hrcast : (r : ℤ) = -h := Int.toNat_of_nonneg (by omega)
    have hrpos : 0 < r := by
      have : (0 : ℤ) < (r : ℤ) := by rw [hrcast]; omega
      exact_mod_cast this
    have hrle : r ≤ d + 1 := by
      rw [← Int.ofNat_le]
      push_cast
      omega
    have hr2 : r ≤ 2 * L := by omega
    have ht : t < 2 * L := by
      dsimp [t]
      omega
    have hn : 2 * P - 1 < 2 * P := by omega
    have hv := multiJumps_block_value (d := d) (n := 2 * P - 1)
      (t := t) (by omega : 0 < L) hn ht
    have hindex :
        ((4 * L * P : ℕ) : ℤ) + h =
          ((2 * L * (2 * P - 1) + t : ℕ) : ℤ) := by
      dsimp [t]
      push_cast
      rw [Nat.cast_sub (by omega : 1 ≤ 2 * P),
        Nat.cast_sub hr2, hrcast]
      push_cast
      ring
    rw [hindex, hv, abs_mul, abs_neg_one_pow, one_mul]
    dsimp [t]
    have hb := abs_blockPoly_right_endpoint_le d (d + 1) L r
      hd heven (by omega) hrle hr2
    unfold realBlock boundaryAmplitudeUpper
    exact_mod_cast hb

theorem left_endpoint_blockAmplitude_le {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hP : 0 < P) (s : ℤ) (hs0 : -(d + 1 : ℕ) ≤ s)
    (hs1 : s ≤ -1) :
    blockAmplitude (d + 2) (multiJumps d L P : Seq) s ≤
      boundaryAmplitudeUpper d L := by
  unfold blockAmplitude
  simp only [dif_pos (by omega : 0 < d + 2)]
  apply Finset.sup'_le
  intro q hq
  simp only [Finset.mem_range] at hq
  apply multiJumps_left_endpoint_value_le hd heven hL hP
  · push_cast
    omega
  · push_cast
    omega

theorem right_endpoint_blockAmplitude_le {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hP : 0 < P) (s : ℤ) (hs0 : -(d + 1 : ℕ) ≤ s)
    (hs1 : s ≤ -1) :
    blockAmplitude (d + 2) (multiJumps d L P : Seq)
      ((4 * L * P : ℕ) + s) ≤ boundaryAmplitudeUpper d L := by
  unfold blockAmplitude
  simp only [dif_pos (by omega : 0 < d + 2)]
  apply Finset.sup'_le
  intro q hq
  simp only [Finset.mem_range] at hq
  rw [show ((4 * L * P : ℕ) : ℤ) + s + (q : ℤ) =
      ((4 * L * P : ℕ) : ℤ) + (s + (q : ℤ)) by ring]
  apply multiJumps_right_endpoint_value_le hd heven hL hP
  · push_cast
    omega
  · push_cast
    omega

theorem left_endpoint_difference_le {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hP : 0 < P) (s : ℤ) (hs0 : -(d + 1 : ℕ) ≤ s)
    (hs1 : s ≤ -1) :
    |cdiffIter (d + 1) (multiJumps d L P) s| ≤
      2 ^ (d + 1) * boundaryAmplitudeUpper d L := by
  rw [cdiffIter_apply]
  apply abs_diffIter_le_of_block
  · unfold boundaryAmplitudeUpper
    exact mul_nonneg (by exact_mod_cast endpointUpper_nonneg d (d + 1))
      (by positivity)
  · intro q hq
    exact multiJumps_left_endpoint_value_le hd heven hL hP (s + q)
      (by push_cast; omega) (by push_cast; omega)

theorem right_endpoint_difference_le {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hP : 0 < P) (s : ℤ) (hs0 : -(d + 1 : ℕ) ≤ s)
    (hs1 : s ≤ -1) :
    |cdiffIter (d + 1) (multiJumps d L P)
        ((4 * L * P : ℕ) + s)| ≤
      2 ^ (d + 1) * boundaryAmplitudeUpper d L := by
  rw [cdiffIter_apply]
  apply abs_diffIter_le_of_block
  · unfold boundaryAmplitudeUpper
    exact mul_nonneg (by exact_mod_cast endpointUpper_nonneg d (d + 1))
      (by positivity)
  · intro q hq
    rw [show ((4 * L * P : ℕ) : ℤ) + s + (q : ℤ) =
      ((4 * L * P : ℕ) : ℤ) + (s + (q : ℤ)) by ring]
    exact multiJumps_right_endpoint_value_le hd heven hL hP (s + q)
      (by push_cast; omega) (by push_cast; omega)

/-! ## Support of the highest difference -/

def realBlockPolynomialEval (d L : ℕ) (x : ℤ) : ℝ :=
  ((blockPolynomial d L).eval (x : ℚ) : ℚ)

theorem realBlock_eq_polynomial (d L t : ℕ) :
    realBlock d L t = realBlockPolynomialEval d L t := by
  unfold realBlockPolynomialEval realBlock
  rw [blockPolynomial_eval]
  norm_num

theorem realBlockPolynomial_high_difference (d L : ℕ) (x : ℤ) :
    diffIter (d + 1) (realBlockPolynomialEval d L) x = 0 := by
  have hbridge :
      diffIter (d + 1) (realBlockPolynomialEval d L) x =
        (((fwdDiff (1 : ℚ))^[d + 1]
          (blockPolynomial d L).eval (x : ℚ) : ℚ) : ℝ) := by
    rw [diffIter_eq_fwdDiff, fwdDiff_iter_eq_sum_shift,
      fwdDiff_iter_eq_sum_shift]
    unfold realBlockPolynomialEval
    push_cast
    apply Finset.sum_congr rfl
    intro q hq
    simp only [zsmul_eq_mul, nsmul_eq_mul, mul_one]
    norm_cast
  rw [hbridge, blockPolynomial_high_difference]
  norm_num

theorem multiJumps_interior_difference_zero {d L P n t : ℕ}
    (hL : 0 < L) (hn : n < 2 * P)
    (ht : t + (d + 1) < 2 * L) :
    cdiffIter (d + 1) (multiJumps d L P)
      (2 * L * n + t : ℕ) = 0 := by
  rw [cdiffIter_apply]
  let B : ℤ := (2 * L * n : ℕ)
  let c : ℝ := (-1 : ℝ) ^ n
  let g : Seq := fun x => c * realBlockPolynomialEval d L (x - B)
  have hcongr :
      diffIter (d + 1) (multiJumps d L P : Seq) (B + t) =
        diffIter (d + 1) g (B + t) := by
    apply diffIter_congr_of_block
    intro q hq
    have htq : t + q < 2 * L := by omega
    have hv := multiJumps_block_value (d := d) (n := n) (t := t + q)
      hL hn htq
    change multiJumps d L P (B + (t : ℤ) + (q : ℤ)) =
      g (B + (t : ℤ) + (q : ℤ))
    rw [show B + (t : ℤ) + (q : ℤ) =
        ((2 * L * n + (t + q) : ℕ) : ℤ) by
          dsimp [B]
          push_cast
          ring,
      hv, realBlock_eq_polynomial]
    simp [g, c, B]
  change diffIter (d + 1) (multiJumps d L P : Seq) (B + t) = 0
  rw [hcongr, diffIter_mul_comp_sub,
    realBlockPolynomial_high_difference, mul_zero]

def boundaryNeighborhood (d L n : ℕ) : Finset ℤ :=
  Finset.Ico (((2 * L * n : ℕ) : ℤ) - (d + 1 : ℕ))
    ((2 * L * n : ℕ) : ℤ)

def boundaryNeighborhoods (d L P : ℕ) : Finset ℤ :=
  (Finset.range (2 * P + 1)).biUnion (boundaryNeighborhood d L)

theorem highestDifference_support_subset {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L) :
    (cdiffIter (d + 1) (multiJumps d L P)).support ⊆
      boundaryNeighborhoods d L P := by
  intro j hj
  have hjne := Finsupp.mem_support_iff.mp hj
  by_cases hjneg : j < 0
  · have hjlow : -((d + 1 : ℕ) : ℤ) ≤ j := by
      by_contra hlow
      have hz :
          cdiffIter (d + 1) (multiJumps d L P) j = 0 := by
        rw [cdiffIter_apply]
        have hc := diffIter_congr_of_block (d + 1)
          (multiJumps d L P : Seq) 0 j (by
            intro q hq
            rw [Pi.zero_apply, multiJumps_zero_of_neg]
            push_cast
            omega)
        rw [hc]
        exact diffIter_zero_seq (d + 1) j
      exact hjne hz
    unfold boundaryNeighborhoods
    rw [Finset.mem_biUnion]
    refine ⟨0, by simp, ?_⟩
    simp only [boundaryNeighborhood, Finset.mem_Ico]
    constructor <;> simp_all
  · have hj0 : 0 ≤ j := le_of_not_gt hjneg
    by_cases hjend : (4 * L * P : ℕ) ≤ j
    · have hz :
          cdiffIter (d + 1) (multiJumps d L P) j = 0 := by
        rw [cdiffIter_apply]
        have hc := diffIter_congr_of_block (d + 1)
          (multiJumps d L P : Seq) 0 j (by
            intro q hq
            rw [Pi.zero_apply, multiJumps_zero_of_ge]
            push_cast at hjend ⊢
            omega)
        rw [hc]
        exact diffIter_zero_seq (d + 1) j
      exact (hjne hz).elim
    · have hjlt : j < (4 * L * P : ℕ) := by omega
      let r := Int.toNat j
      let n := r / (2 * L)
      let t := r % (2 * L)
      have hrcast : (r : ℤ) = j := Int.toNat_of_nonneg hj0
      have hden : 0 < 2 * L := by omega
      have ht : t < 2 * L := Nat.mod_lt _ hden
      have hdecomp : 2 * L * n + t = r := by
        simpa [n, t, mul_comm] using (Nat.div_add_mod r (2 * L))
      have hn : n < 2 * P := by
        rw [Nat.div_lt_iff_lt_mul hden]
        have : r < 4 * L * P := by
          rw [← hrcast] at hjlt
          exact_mod_cast hjlt
        nlinarith
      by_cases hinterior : t + (d + 1) < 2 * L
      · have hz := multiJumps_interior_difference_zero
          (d := d) (L := L) (P := P) (by omega) hn hinterior
        have hindex : (2 * L * n + t : ℕ) = j := by
          rw [hdecomp]
          exact_mod_cast hrcast
        rw [hindex] at hz
        exact (hjne hz).elim
      · unfold boundaryNeighborhoods
        rw [Finset.mem_biUnion]
        refine ⟨n + 1, by
          simp only [Finset.mem_range]
          omega, ?_⟩
        simp only [boundaryNeighborhood, Finset.mem_Ico]
        have hindex : ((2 * L * n + t : ℕ) : ℤ) = j := by
          rw [hdecomp]
          exact hrcast
        constructor
        · rw [← hindex]
          have hnot : 2 * L ≤ t + (d + 1) := by omega
          push_cast
          push_cast at hnot
          nlinarith
        · rw [← hindex]
          push_cast
          have htZ : (t : ℤ) < 2 * (L : ℤ) := by exact_mod_cast ht
          calc
            2 * (L : ℤ) * (n : ℤ) + (t : ℤ) <
                2 * (L : ℤ) * (n : ℤ) + 2 * (L : ℤ) :=
              by simpa [add_comm] using
                add_lt_add_left htZ (2 * (L : ℤ) * (n : ℤ))
            _ = 2 * (L : ℤ) * ((n : ℤ) + 1) := by ring

theorem boundaryNeighborhood_pairwiseDisjoint {d L P : ℕ}
    (hL : d + 1 ≤ L) :
    ((Finset.range (2 * P + 1) : Set ℕ)).PairwiseDisjoint
      (boundaryNeighborhood d L) := by
  intro n hn m hm hnm
  change Disjoint (boundaryNeighborhood d L n)
    (boundaryNeighborhood d L m)
  rw [Finset.disjoint_left]
  intro j hjn hjm
  simp only [boundaryNeighborhood, Finset.mem_Ico] at hjn hjm
  by_cases hnm' : n < m
  · have hstep : n + 1 ≤ m := by omega
    have hnat : 2 * L * (n + 1) ≤ 2 * L * m :=
      Nat.mul_le_mul_left (2 * L) hstep
    have hcast : ((2 * L * (n + 1) : ℕ) : ℤ) ≤
        ((2 * L * m : ℕ) : ℤ) := by exact_mod_cast hnat
    push_cast at hcast hjn hjm
    nlinarith
  · have hstep : m + 1 ≤ n := by omega
    have hnat : 2 * L * (m + 1) ≤ 2 * L * n :=
      Nat.mul_le_mul_left (2 * L) hstep
    have hcast : ((2 * L * (m + 1) : ℕ) : ℤ) ≤
        ((2 * L * n : ℕ) : ℤ) := by exact_mod_cast hnat
    push_cast at hcast hjn hjm
    nlinarith

def boundaryEnergy (d L P n : ℕ) : ℝ :=
  ∑ j ∈ boundaryNeighborhood d L n,
    blockAmplitude (d + 2) (multiJumps d L P : Seq) j *
      |cdiffIter (d + 1) (multiJumps d L P) j|

theorem boundaryAmplitudeUpper_nonneg (d L : ℕ) :
    0 ≤ boundaryAmplitudeUpper d L := by
  unfold boundaryAmplitudeUpper
  exact mul_nonneg (by exact_mod_cast endpointUpper_nonneg d (d + 1))
    (by positivity)

theorem boundaryNeighborhood_card (d L n : ℕ) :
    (boundaryNeighborhood d L n).card = d + 1 := by
  simp [boundaryNeighborhood, Int.card_Ico]

theorem internal_boundaryEnergy_le {d L P n : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hn0 : 1 ≤ n) (hnP : n < 2 * P) :
    boundaryEnergy d L P n ≤
      (d + 1) *
        (boundaryAmplitudeUpper d L * (junctionUpper d : ℝ)) := by
  unfold boundaryEnergy
  calc
    (∑ j ∈ boundaryNeighborhood d L n,
        blockAmplitude (d + 2) (multiJumps d L P : Seq) j *
          |cdiffIter (d + 1) (multiJumps d L P) j|) ≤
      ∑ _j ∈ boundaryNeighborhood d L n,
        boundaryAmplitudeUpper d L * (junctionUpper d : ℝ) := by
          apply Finset.sum_le_sum
          intro j hj
          simp only [boundaryNeighborhood, Finset.mem_Ico] at hj
          push_cast at hj
          let s := j - (2 * L * n : ℕ)
          have hs0 : -((d + 1 : ℕ) : ℤ) ≤ s := by
            dsimp [s]
            push_cast
            nlinarith
          have hs1 : s ≤ -1 := by
            dsimp [s]
            push_cast
            nlinarith
          have hA := internal_blockAmplitude_le hd heven hL hn0 hnP s hs0 hs1
          have hD := multiJumps_internal_difference_le hd heven hL hn0 hnP s hs0 hs1
          have hjform : j = ((2 * L * n : ℕ) : ℤ) + s := by
            dsimp [s]
            ring
          rw [hjform]
          exact mul_le_mul hA hD (abs_nonneg _)
            (boundaryAmplitudeUpper_nonneg d L)
    _ = (d + 1) *
        (boundaryAmplitudeUpper d L * (junctionUpper d : ℝ)) := by
      rw [Finset.sum_const, nsmul_eq_mul, boundaryNeighborhood_card]
      norm_cast

def endpointEnergyUpper (d L : ℕ) : ℝ :=
  boundaryAmplitudeUpper d L *
    (2 ^ (d + 1) * boundaryAmplitudeUpper d L)

theorem left_boundaryEnergy_le {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hP : 0 < P) :
    boundaryEnergy d L P 0 ≤ (d + 1) * endpointEnergyUpper d L := by
  unfold boundaryEnergy
  calc
    (∑ j ∈ boundaryNeighborhood d L 0,
        blockAmplitude (d + 2) (multiJumps d L P : Seq) j *
          |cdiffIter (d + 1) (multiJumps d L P) j|) ≤
      ∑ _j ∈ boundaryNeighborhood d L 0, endpointEnergyUpper d L := by
        apply Finset.sum_le_sum
        intro j hj
        simp only [boundaryNeighborhood, Finset.mem_Ico] at hj
        have hA := left_endpoint_blockAmplitude_le hd heven hL hP j
          (by simpa using hj.1) (by omega)
        have hD := left_endpoint_difference_le hd heven hL hP j
          (by simpa using hj.1) (by omega)
        unfold endpointEnergyUpper
        exact mul_le_mul hA hD (abs_nonneg _)
          (boundaryAmplitudeUpper_nonneg d L)
    _ = (d + 1) * endpointEnergyUpper d L := by
      rw [Finset.sum_const, nsmul_eq_mul, boundaryNeighborhood_card]
      norm_cast

theorem right_boundaryEnergy_le {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hP : 0 < P) :
    boundaryEnergy d L P (2 * P) ≤
      (d + 1) * endpointEnergyUpper d L := by
  unfold boundaryEnergy
  calc
    (∑ j ∈ boundaryNeighborhood d L (2 * P),
        blockAmplitude (d + 2) (multiJumps d L P : Seq) j *
          |cdiffIter (d + 1) (multiJumps d L P) j|) ≤
      ∑ _j ∈ boundaryNeighborhood d L (2 * P),
        endpointEnergyUpper d L := by
        apply Finset.sum_le_sum
        intro j hj
        simp only [boundaryNeighborhood, Finset.mem_Ico] at hj
        push_cast at hj
        let s := j - (4 * L * P : ℕ)
        have hs0 : -((d + 1 : ℕ) : ℤ) ≤ s := by
          dsimp [s]
          push_cast
          nlinarith
        have hs1 : s ≤ -1 := by
          dsimp [s]
          push_cast
          nlinarith
        have hA := right_endpoint_blockAmplitude_le hd heven hL hP s hs0 hs1
        have hD := right_endpoint_difference_le hd heven hL hP s hs0 hs1
        have hjform : j = ((4 * L * P : ℕ) : ℤ) + s := by
          dsimp [s]
          ring
        rw [hjform]
        exact mul_le_mul hA hD (abs_nonneg _)
          (boundaryAmplitudeUpper_nonneg d L)
    _ = (d + 1) * endpointEnergyUpper d L := by
      rw [Finset.sum_const, nsmul_eq_mul, boundaryNeighborhood_card]
      norm_cast

theorem cEnergy_multiJumps_le_boundary_sum {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L) :
    cEnergy (d + 2) (multiJumps d L P) ≤
      ∑ n ∈ Finset.range (2 * P + 1), boundaryEnergy d L P n := by
  unfold cEnergy
  change (∑ j ∈ (cdiffIter (d + 1) (multiJumps d L P)).support,
      blockAmplitude (d + 2) (multiJumps d L P : Seq) j *
        |cdiffIter (d + 1) (multiJumps d L P) j|) ≤ _
  calc
    _ ≤ ∑ j ∈ (Finset.range (2 * P + 1)).biUnion
          (boundaryNeighborhood d L),
        blockAmplitude (d + 2) (multiJumps d L P : Seq) j *
          |cdiffIter (d + 1) (multiJumps d L P) j| := by
      apply Finset.sum_le_sum_of_subset_of_nonneg
        (highestDifference_support_subset hd heven hL)
      intro j hj hnot
      exact mul_nonneg (blockAmplitude_nonneg _ _ _)
        (abs_nonneg _)
    _ = ∑ n ∈ Finset.range (2 * P + 1),
          ∑ j ∈ boundaryNeighborhood d L n,
            blockAmplitude (d + 2) (multiJumps d L P : Seq) j *
              |cdiffIter (d + 1) (multiJumps d L P) j| := by
      rw [Finset.sum_biUnion (boundaryNeighborhood_pairwiseDisjoint hL)]
    _ = ∑ n ∈ Finset.range (2 * P + 1),
          boundaryEnergy d L P n := by rfl

theorem cEnergy_multiJumps_upper {d L P : ℕ}
    (hd : 0 < d) (heven : Even d) (hL : d + 1 ≤ L)
    (hP : 0 < P) :
    cEnergy (d + 2) (multiJumps d L P) ≤
      ((2 * P - 1 : ℕ) : ℝ) * (d + 1) *
          (boundaryAmplitudeUpper d L * (junctionUpper d : ℝ)) +
        2 * (d + 1) * endpointEnergyUpper d L := by
  let f := boundaryEnergy d L P
  have hbase := cEnergy_multiJumps_le_boundary_sum
    (d := d) (L := L) (P := P) hd heven hL
  have htwo : 2 * P = (2 * P - 1) + 1 := by omega
  have hsplit :
      (∑ n ∈ Finset.range (2 * P + 1), f n) =
        (∑ q ∈ Finset.range (2 * P - 1), f (q + 1)) +
          f (2 * P) + f 0 := by
    calc
      (∑ n ∈ Finset.range (2 * P + 1), f n) =
          (∑ q ∈ Finset.range (2 * P), f (q + 1)) + f 0 :=
        Finset.sum_range_succ' f (2 * P)
      _ = ((∑ q ∈ Finset.range (2 * P - 1), f (q + 1)) +
            f ((2 * P - 1) + 1)) + f 0 := by
        conv_lhs =>
          lhs
          rw [htwo, Finset.sum_range_succ]
      _ = (∑ q ∈ Finset.range (2 * P - 1), f (q + 1)) +
            f (2 * P) + f 0 := by rw [← htwo]
  have hinter :
      (∑ q ∈ Finset.range (2 * P - 1), f (q + 1)) ≤
        ((2 * P - 1 : ℕ) : ℝ) * ((d + 1) *
          (boundaryAmplitudeUpper d L * (junctionUpper d : ℝ))) := by
    calc
      (∑ q ∈ Finset.range (2 * P - 1), f (q + 1)) ≤
          ∑ _q ∈ Finset.range (2 * P - 1),
            (d + 1) *
              (boundaryAmplitudeUpper d L * (junctionUpper d : ℝ)) := by
        apply Finset.sum_le_sum
        intro q hq
        simp only [Finset.mem_range] at hq
        exact internal_boundaryEnergy_le hd heven hL (by omega) (by omega)
      _ = ((2 * P - 1 : ℕ) : ℝ) * ((d + 1) *
          (boundaryAmplitudeUpper d L * (junctionUpper d : ℝ))) := by
        rw [Finset.sum_const, nsmul_eq_mul]
        rw [Finset.card_range]
  have hright := right_boundaryEnergy_le hd heven hL hP
  have hleft := left_boundaryEnergy_le hd heven hL hP
  change f (2 * P) ≤ (d + 1) * endpointEnergyUpper d L at hright
  change f 0 ≤ (d + 1) * endpointEnergyUpper d L at hleft
  calc
    cEnergy (d + 2) (multiJumps d L P) ≤
        ∑ n ∈ Finset.range (2 * P + 1), f n := hbase
    _ = (∑ q ∈ Finset.range (2 * P - 1), f (q + 1)) +
          f (2 * P) + f 0 := hsplit
    _ ≤ ((2 * P - 1 : ℕ) : ℝ) * ((d + 1) *
          (boundaryAmplitudeUpper d L * (junctionUpper d : ℝ))) +
        (d + 1) * endpointEnergyUpper d L +
        (d + 1) * endpointEnergyUpper d L := by gcongr
    _ = ((2 * P - 1 : ℕ) : ℝ) * (d + 1) *
          (boundaryAmplitudeUpper d L * (junctionUpper d : ℝ)) +
        2 * (d + 1) * endpointEnergyUpper d L := by ring

end ENOTV