CoupledNormalization.lean

1 import CoupledEvaluation
2 import Mathlib.Data.Nat.Choose.Basic
3 
4 /-!
5 # Factorial normalization of the coupled derivative recurrence
6 
7 The mathematical recurrence underlying `coupled_odd_derivatives` is defined
8 over exact rationals. Its order-(2n+1) arrays, divided by (2n+1)!, equal the
9 ordinary Taylor recurrence at every index, for arbitrary rational scales.
10 Finite inclusive prefixes therefore yield the certified polynomials.
11 
12 This is a mathematical recurrence and finite-array specification. It does
13 not certify a Python interpreter, Fraction, mutation, guards or exceptions.
14 -/
15 
16 set_option autoImplicit false
17 set_option maxHeartbeats 2500000
18 
19 noncomputable section
20 
21 namespace Komorebi.SpecialFunctions.Coupled.Scales
22 
23 def oddFactorial (n : ℕ) : ℚ := ((2*n+1).factorial : ℚ)
24 
25 theorem oddFactorial_ne_zero (n : ℕ) : oddFactorial n ≠ 0 := by
26   unfold oddFactorial
27   exact_mod_cast Nat.factorial_ne_zero (2*n+1)
28 
29 theorem oddFactorial_zero : oddFactorial 0 = 1 := by
30   norm_num [oddFactorial]
31 
32 theorem oddFactorial_succ (n : ℕ) :
33     oddFactorial (n+1) = (2*(n : ℚ)+3) * ((2*n+2).factorial : ℚ) := by
34   have he : 2*(n+1)+1 = (2*n+2)+1 := by omega
35   rw [oddFactorial, he, Nat.factorial_succ]
36   push_cast
37   ring
38 
39 theorem binomial_factorial_product (n i : ℕ) (hi : i ≤ n) :
40     ((2*n+2).choose (2*i+1) : ℚ) * oddFactorial i * oddFactorial (n-i) =
41       ((2*n+2).factorial : ℚ) := by
42   have h := Nat.choose_mul_factorial_mul_factorial (show 2*i+1 ≤ 2*n+2 by omega)
43   have he : 2*n+2-(2*i+1) = 2*(n-i)+1 := by omega
44   rw [he] at h
45   unfold oddFactorial
46   exact_mod_cast h
47 
48 theorem binomial_factorial_quotient (n i : ℕ) (hi : i ≤ n) :
49     ((2*n+2).choose (2*i+1) : ℚ) / oddFactorial (n+1) =
50       1 / ((2*(n : ℚ)+3) * oddFactorial i * oddFactorial (n-i)) := by
51   have h := binomial_factorial_product n i hi
52   have hf := oddFactorial_succ n
53   have hn : (2*(n : ℚ)+3) ≠ 0 := by positivity
54   apply (div_eq_div_iff (oddFactorial_ne_zero (n+1))
55     (mul_ne_zero (mul_ne_zero hn (oddFactorial_ne_zero i)) (oddFactorial_ne_zero (n-i)))).mpr
56   calc
57     _ = (2*(n : ℚ)+3) *
58         (((2*n+2).choose (2*i+1) : ℚ) * oddFactorial i * oddFactorial (n-i)) := by ring
59     _ = _ := by rw [h, hf]; ring
60 
61 def derivativeRow (q r : ℚ) (c : ℕ → ℚ × ℚ) (n : ℕ) : ℚ :=
62   ∑ i : Fin (n+1), ((2*n+2).choose (2*i.val+1) : ℚ) *
63     (c i.val).1 * (c (n-i.val)).2 * q^(2*i.val+1) * r^(2*(n-i.val)+1)
64 
65 def RationalScales.derivativeNext (p : RationalScales) (c : ℕ → ℚ × ℚ) (n : ℕ) : ℚ × ℚ :=
66   (derivativeRow p.q p.r c n, derivativeRow p.a p.b c n)
67 
68 def RationalScales.derivatives (p : RationalScales) : ℕ → ℚ × ℚ
69   | 0 => (1, 1)
70   | n+1 =>
71     ((∑ i : Fin (n+1), ((2*n+2).choose (2*i.val+1) : ℚ) *
72       (p.derivatives i.val).1 * (p.derivatives (n-i.val)).2 *
73       p.q^(2*i.val+1) * p.r^(2*(n-i.val)+1)),
74      (∑ i : Fin (n+1), ((2*n+2).choose (2*i.val+1) : ℚ) *
75       (p.derivatives i.val).1 * (p.derivatives (n-i.val)).2 *
76       p.a^(2*i.val+1) * p.b^(2*(n-i.val)+1)))
77 termination_by n => n
78 decreasing_by all_goals omega
79 
80 theorem rational_derivatives_zero (p : RationalScales) : p.derivatives 0 = (1,1) := by
81   rw [RationalScales.derivatives.eq_def]
82 
83 theorem rational_derivatives_succ (p : RationalScales) (n : ℕ) :
84     p.derivatives (n+1) = p.derivativeNext p.derivatives n := by
85   rw [RationalScales.derivatives.eq_def]
86   rfl
87 
88 def RationalScales.IsDerivativeSolution (p : RationalScales) (c : ℕ → ℚ × ℚ) : Prop :=
89   c 0 = (1,1) ∧ ∀ n, c (n+1) = p.derivativeNext c n
90 
91 theorem rational_derivatives_solution (p : RationalScales) :
92     p.IsDerivativeSolution p.derivatives :=
93   ⟨rational_derivatives_zero p, rational_derivatives_succ p⟩
94 
95 theorem derivative_solution_unique (p : RationalScales) {c d : ℕ → ℚ × ℚ}
96     (hc : p.IsDerivativeSolution c) (hd : p.IsDerivativeSolution d) : c = d := by
97   funext n
98   induction n using Nat.strong_induction_on with
99   | h n ih =>
100     cases n with
101     | zero => exact hc.1.trans hd.1.symm
102     | succ n =>
103       rw [hc.2 n, hd.2 n]
104       apply Prod.ext
105       all_goals
106         simp only [RationalScales.derivativeNext, derivativeRow]
107         apply Finset.sum_congr rfl
108         intro i _
109         have hi : i.val < n+1 := i.isLt
110         have hj : n-i.val < n+1 := by omega
111         simp only [ih i.val hi, ih (n-i.val) hj]
112 
113 theorem derivativeRow_factorialScale (q r : ℚ) (c : ℕ → ℚ × ℚ) (n : ℕ) :
114     derivativeRow q r (fun k => (oddFactorial k * (c k).1, oddFactorial k * (c k).2)) n =
115       ((2*n+2).factorial : ℚ) *
116         ∑ i : Fin (n+1), q^(2*i.val+1) * r^(2*(n-i.val)+1) * (c i.val).1 * (c (n-i.val)).2 := by
117   rw [derivativeRow, Finset.mul_sum]
118   apply Finset.sum_congr rfl
119   intro i _
120   have hi : i.val ≤ n := by omega
121   have h := binomial_factorial_product n i.val hi
122   dsimp only
123   calc
124     _ = (((2*n+2).choose (2*i.val+1) : ℚ) * oddFactorial i.val * oddFactorial (n-i.val)) *
125         (q^(2*i.val+1) * r^(2*(n-i.val)+1) * (c i.val).1 * (c (n-i.val)).2) := by ring
126     _ = _ := by rw [h]
127 
128 theorem factorial_scaled_coefficients_solution (p : RationalScales) :
129     p.IsDerivativeSolution
130       (fun n => (oddFactorial n * (p.coefficients n).1, oddFactorial n * (p.coefficients n).2)) := by
131   constructor
132   · dsimp only
133     have hz : p.coefficients 0 = (1,1) := by rw [RationalScales.coefficients.eq_def]
134     rw [hz, oddFactorial_zero]
135     norm_num
136   · intro n
137     dsimp only
138     have hn : (2*(n : ℚ)+3) ≠ 0 := by positivity
139     have hc := RationalScales.coefficients.eq_def p (n+1)
140     rw [hc]
141     apply Prod.ext
142     all_goals
143       simp only [RationalScales.derivativeNext, derivativeRow_factorialScale,
144         oddFactorial_succ]
145       field_simp [hn]
146 
147 theorem rational_derivatives_factorial (p : RationalScales) (n : ℕ) :
148     p.derivatives n =
149       (oddFactorial n * (p.coefficients n).1, oddFactorial n * (p.coefficients n).2) :=
150   congrFun (derivative_solution_unique p (rational_derivatives_solution p)
151     (factorial_scaled_coefficients_solution p)) n
152 
153 theorem derivativeNext_congr (p : RationalScales) {c d : ℕ → ℚ × ℚ} (n : ℕ)
154     (h : ∀ j, j ≤ n → c j = d j) : p.derivativeNext c n = p.derivativeNext d n := by
155   apply Prod.ext
156   all_goals
157     simp only [RationalScales.derivativeNext, derivativeRow]
158     apply Finset.sum_congr rfl
159     intro i _
160     rw [h i.val (by omega), h (n-i.val) (Nat.sub_le _ _)]
161 
162 /-- Paired-array append model; both rows use only entries already stored. -/
163 def RationalScales.derivativeArray (p : RationalScales) : ℕ → List (ℚ × ℚ)
164   | 0 => [(1,1)]
165   | n+1 =>
166     let prior := p.derivativeArray n
167     prior ++ [p.derivativeNext (fun k => (prior[k]?).getD (0,0)) n]
168 
169 theorem derivative_array_eq (p : RationalScales) (last : ℕ) :
170     p.derivativeArray last = List.ofFn (fun i : Fin (last+1) => p.derivatives i.val) := by
171   induction last with
172   | zero => simp [RationalScales.derivativeArray, List.ofFn_succ, rational_derivatives_zero]
173   | succ n ih =>
174     rw [RationalScales.derivativeArray, ih]
175     have hnext : p.derivativeNext
176         (fun k => ((List.ofFn (fun i : Fin (n+1) => p.derivatives i.val))[k]?).getD (0,0)) n =
177         p.derivatives (n+1) := by
178       rw [rational_derivatives_succ]
179       apply derivativeNext_congr
180       intro j hj
181       simp only [List.getElem?_ofFn, dite_eq_left (show j < n+1 by omega), Option.getD_some]
182     rw [hnext]
183     exact (List.ofFn_succ_last (f := fun j : Fin (n+1+1) => p.derivatives j.val)).symm
184 
185 theorem derivative_array_length (p : RationalScales) (last : ℕ) :
186     (p.derivativeArray last).length = last+1 := by
187   rw [derivative_array_eq, List.length_ofFn]
188 
189 theorem derivative_array_entry (p : RationalScales) (last n : ℕ) (hn : n ≤ last) :
190     (p.derivativeArray last)[n]? = some (p.derivatives n) := by
191   rw [derivative_array_eq]
192   simp only [List.getElem?_ofFn, dite_eq_left (show n < last+1 by omega)]
193 
194 def RationalScales.derivativeComponent (p : RationalScales) (i : Fin 2) (n : ℕ) : ℚ :=
195   ![(p.derivatives n).1, (p.derivatives n).2] i
196 
197 theorem derivative_component_factorial (p : RationalScales) (i : Fin 2) (n : ℕ) :
198     p.derivativeComponent i n = oddFactorial n * p.componentCoefficient i n := by
199   fin_cases i <;> simp [RationalScales.derivativeComponent, RationalScales.componentCoefficient,
200     rational_derivatives_factorial]
201 
202 theorem derivative_component_normalized (p : RationalScales) (i : Fin 2) (n : ℕ) :
203     p.derivativeComponent i n / oddFactorial n = p.componentCoefficient i n := by
204   rw [derivative_component_factorial]
205   exact mul_div_cancel_left₀ _ (oddFactorial_ne_zero n)
206 
207 theorem derivative_normalized_cast (p : RationalScales) (i : Fin 2) (n : ℕ) :
208     ((p.derivativeComponent i n / oddFactorial n : ℚ) : ℝ) =
209       p.toReal.componentCoefficient i n := by
210   rw [derivative_component_normalized, rational_component_coefficient_cast]
211 
212 def RationalScales.ordinaryArray (p : RationalScales) (i : Fin 2) (last : ℕ) : List ℚ :=
213   (p.derivativeArray last).mapIdx (fun n c => ![c.1, c.2] i / oddFactorial n)
214 
215 theorem ordinary_array_eq (p : RationalScales) (i : Fin 2) (last : ℕ) :
216     p.ordinaryArray i last = List.ofFn (fun n : Fin (last+1) => p.componentCoefficient i n.val) := by
217   rw [RationalScales.ordinaryArray, derivative_array_eq]
218   apply List.ext_getElem
219   · simp
220   · intro n hn hm
221     simp only [List.getElem_mapIdx, List.getElem_ofFn]
222     exact derivative_component_normalized p i n
223 
224 theorem ordinary_array_length (p : RationalScales) (i : Fin 2) (last : ℕ) :
225     (p.ordinaryArray i last).length = last+1 := by
226   rw [ordinary_array_eq, List.length_ofFn]
227 
228 theorem ordinary_array_entry (p : RationalScales) (i : Fin 2) (last n : ℕ) (hn : n ≤ last) :
229     (p.ordinaryArray i last)[n]? = some (p.componentCoefficient i n) := by
230   rw [ordinary_array_eq]
231   simp only [List.getElem?_ofFn, dite_eq_left (show n < last+1 by omega)]
232 
233 def RationalScales.derivativePolynomial (p : RationalScales) (i : Fin 2) (N : ℕ) (x : ℚ) : ℚ :=
234   ∑ n ∈ Finset.range (N+1), (p.derivativeComponent i n / oddFactorial n) * x^(2*n+1)
235 
236 theorem derivative_polynomial_eq (p : RationalScales) (i : Fin 2) (N : ℕ) (x : ℚ) :
237     p.derivativePolynomial i N x = p.polynomial i N x := by
238   simp only [RationalScales.derivativePolynomial, RationalScales.polynomial,
239     derivative_component_normalized]
240 
241 def RationalScales.arrayPolynomial (p : RationalScales) (i : Fin 2) (N : ℕ) (x : ℚ) : ℚ :=
242   ((p.ordinaryArray i N).mapIdx (fun n c => c*x^(2*n+1))).sum
243 
244 theorem array_polynomial_eq (p : RationalScales) (i : Fin 2) (N : ℕ) (x : ℚ) :
245     p.arrayPolynomial i N x = p.polynomial i N x := by
246   have hlist : (p.ordinaryArray i N).mapIdx (fun n c => c*x^(2*n+1)) =
247       List.ofFn (fun n : Fin (N+1) => p.componentCoefficient i n.val * x^(2*n.val+1)) := by
248     rw [ordinary_array_eq]
249     apply List.ext_getElem
250     · simp
251     · intro n hn hm
252       simp only [List.getElem_mapIdx, List.getElem_ofFn]
253   rw [RationalScales.arrayPolynomial, hlist, Fin.sum_ofFn, RationalScales.polynomial,
254     Finset.sum_range]
255 
256 theorem array_polynomial_cast (p : RationalScales) (i : Fin 2) (N : ℕ) (x : ℚ) :
257     (p.arrayPolynomial i N x : ℝ) = p.toReal.realOddPolynomial i N (x : ℝ) := by
258   rw [array_polynomial_eq, rational_polynomial_cast]
259 
260 theorem derivative_polynomial_cast (p : RationalScales) (i : Fin 2) (N : ℕ) (x : ℚ) :
261     (p.derivativePolynomial i N x : ℝ) = p.toReal.realOddPolynomial i N (x : ℝ) := by
262   rw [derivative_polynomial_eq, rational_polynomial_cast]
263 
264 theorem derivative_endpoint_enclosure (p : RationalScales) (hs : p.toReal.positive)
265     (hc : p.toReal.q < 1 ∧ p.toReal.r < 1 ∧ p.toReal.a < 1 ∧ p.toReal.b < 1)
266     (x S : ℚ) (hx : 0 ≤ x) (hxS : x < S) (fuel : ℕ) (ht : p.maxScale^fuel*S < 1)
267     (N : ℕ) (i : Fin 2) :
268     (p.derivativePolynomial i N x : ℝ) ≤ p.toReal.realComponent i (x : ℝ) ∧
269       p.toReal.realComponent i (x : ℝ) ≤
270         ((p.derivativePolynomial i N x + rationalDiskMajorant p.maxScale fuel S *
271           (x/S)^(2*N+3) : ℚ) : ℝ) := by
272   rw [derivative_polynomial_eq]
273   exact rational_endpoint_enclosure p hs hc x S hx hxS fuel ht N i
274 
275 theorem array_endpoint_enclosure (p : RationalScales) (hs : p.toReal.positive)
276     (hc : p.toReal.q < 1 ∧ p.toReal.r < 1 ∧ p.toReal.a < 1 ∧ p.toReal.b < 1)
277     (x S : ℚ) (hx : 0 ≤ x) (hxS : x < S) (fuel : ℕ) (ht : p.maxScale^fuel*S < 1)
278     (N : ℕ) (i : Fin 2) :
279     (p.arrayPolynomial i N x : ℝ) ≤ p.toReal.realComponent i (x : ℝ) ∧
280       p.toReal.realComponent i (x : ℝ) ≤
281         ((p.arrayPolynomial i N x + rationalDiskMajorant p.maxScale fuel S *
282           (x/S)^(2*N+3) : ℚ) : ℝ) := by
283   rw [array_polynomial_eq]
284   exact rational_endpoint_enclosure p hs hc x S hx hxS fuel ht N i
285 
286 theorem rational_derivatives_one (p : RationalScales) :
287     p.derivatives 1 = (2*p.q*p.r, 2*p.a*p.b) := by
288   rw [rational_derivatives_succ]
289   simp [RationalScales.derivativeNext, derivativeRow, rational_derivatives_zero]
290 
291 theorem derivative_normalization_first (p : RationalScales) :
292     p.derivativeComponent 0 1 / oddFactorial 1 = p.q*p.r/3 ∧
293       p.derivativeComponent 1 1 / oddFactorial 1 = p.a*p.b/3 := by
294   simp [RationalScales.derivativeComponent, rational_derivatives_one, oddFactorial]
295   constructor <;> ring
296 
297 end Komorebi.SpecialFunctions.Coupled.Scales

Back to source controls