Require Import Reals Omega Psatz.
Require Import mathcomp.ssreflect.ssreflect.
Require Import Rcomplements Hierarchy Continuity Derive.

Definition differentiable_pt_lim (f : R → R → R) (x y : R) (lx ly : R) :=
  ∀ eps : posreal, locally_2d (fun u v ⇒
    Rabs (f u v - f x y - (lx × (u - x) + ly × (v - y))) ≤ eps × Rmax (Rabs (u - x)) (Rabs (v - y))) x y.

Lemma filterdiff_differentiable_pt_lim (f : R → R → R) (x y lx ly : R) :
  filterdiff (fun u : R × R ⇒ f (fst u) (snd u)) (locally (x,y)) (fun u : R × R ⇒ fst u × lx + snd u × ly)
  ↔ differentiable_pt_lim f x y lx ly.

Lemma differentiable_pt_lim_ext : ∀ f1 f2 x y lx ly,
  locally_2d (fun u v ⇒ f1 u v = f2 u v) x y →
  differentiable_pt_lim f1 x y lx ly → differentiable_pt_lim f2 x y lx ly.

Definition differentiable_pt (f : R → R → R) (x y : R) :=
  ∃ lx, ∃ ly, differentiable_pt_lim f x y lx ly.

Lemma differentiable_continuity_pt : ∀ f x y,
  differentiable_pt f x y → continuity_2d_pt f x y.

Lemma differentiable_pt_lim_proj1_0 (f : R → R) (x y l : R) :
  derivable_pt_lim f x l → differentiable_pt_lim (fun u v ⇒ f u) x y l 0.

Lemma differentiable_pt_lim_proj1_1 (f : R → R) (x y l : R) :
  differentiable_pt_lim (fun u v ⇒ f u) x y l 0 → derivable_pt_lim f x l.

Lemma differentiable_pt_lim_unique (f : R → R → R) (x y : R) (lx ly : R) :
  differentiable_pt_lim f x y lx ly
    → Derive (fun x ⇒ f x y) x = lx ∧ Derive (fun y ⇒ f x y) y = ly.

Lemma differentiable_pt_lim_comp : ∀ f1 f2 f3 x y l1x l1y l2x l2y l3x l3y,
  differentiable_pt_lim f1 (f2 x y) (f3 x y) l1x l1y →
  differentiable_pt_lim f2 x y l2x l2y → differentiable_pt_lim f3 x y l3x l3y →
  differentiable_pt_lim (fun u v ⇒ f1 (f2 u v) (f3 u v)) x y
    (l1x × l2x + l1y × l3x) (l1x × l2y + l1y × l3y).

Lemma derivable_pt_lim_comp_2d : ∀ f1 f2 f3 x l1x l1y l2 l3,
  differentiable_pt_lim f1 (f2 x) (f3 x) l1x l1y →
  derivable_pt_lim f2 x l2 → derivable_pt_lim f3 x l3 →
  derivable_pt_lim (fun t ⇒ f1 (f2 t) (f3 t)) x (l1x × l2 + l1y × l3).

Definition partial_derive (m k : nat) (f : R → R → R) : R → R → R :=
  fun x y ⇒ Derive_n (fun t ⇒ Derive_n (fun z ⇒ f t z) k y) m x.

Definition differential (p : nat) (f : R → R → R) (x y dx dy : R) : R :=
  sum_f_R0
    (fun m ⇒
      C p m ×
      partial_derive m (p - m)%nat f x y ×
      dx ^ m × dy ^ (p - m)%nat)
    p.

Definition DL_pol (n : nat) (f : R → R → R) (x y dx dy : R) : R :=
  sum_f_R0
    (fun p ⇒
      differential p f x y dx dy / INR (fact p))
    n.

Lemma partial_derive_ext_loc :
  ∀ f g p q x y,
  locally_2d (fun u v ⇒ f u v = g u v) x y →
  partial_derive p q f x y = partial_derive p q g x y.

Lemma Schwarz_aux :
  ∀ f x y (eps : posreal),
  ( ∀ u v, Rabs (u - x) < eps → Rabs (v - y) < eps →
    ex_derive (fun z : R ⇒ f z v) u ∧
    ex_derive (fun z : R ⇒ Derive (fun t ⇒ f t z) u) v ) →
  ∀ h k, Rabs h < eps → Rabs k < eps →
  let phi k x := f x (y + k) - f x y in
  ∃ u, ∃ v,
  Rabs (u - x) ≤ Rabs h ∧ Rabs (v - y) ≤ Rabs k ∧
  phi k (x + h) - phi k x = h × k × (Derive (fun z ⇒ Derive (fun t ⇒ f t z) u) v).

Lemma Schwarz :
  ∀ (f : R → R → R) x y,
  locally_2d (fun u v ⇒
    ex_derive (fun z : R ⇒ f z v) u ∧
    ex_derive (fun z : R ⇒ f u z) v ∧
    ex_derive (fun z : R ⇒ Derive (fun t ⇒ f z t) v) u ∧
    ex_derive (fun z : R ⇒ Derive (fun t ⇒ f t z) u) v) x y →
  continuity_2d_pt (fun u v ⇒ Derive (fun z ⇒ Derive (fun t ⇒ f z t) v) u) x y →
  continuity_2d_pt (fun u v ⇒ Derive (fun z ⇒ Derive (fun t ⇒ f t z) u) v) x y →
  Derive (fun z ⇒ Derive (fun t ⇒ f z t) y) x = Derive (fun z ⇒ Derive (fun t ⇒ f t z) x) y.

Lemma partial_derive_add_zero: ∀ f p q r s x y,
  (q=0)%nat ∨ (r=0)%nat →
  partial_derive p q (partial_derive r s f) x y
   = partial_derive (p+r) (q+s) f x y.

Fixpoint ex_diff_n f n x y :=
  continuity_2d_pt f x y ∧
  match n with
  | O ⇒ True
  | S n ⇒
    ex_derive (fun z ⇒ f z y) x ∧
    ex_derive (fun z ⇒ f x z) y ∧
    ex_diff_n (fun u v ⇒ Derive (fun z ⇒ f z v) u) n x y ∧
    ex_diff_n (fun u v ⇒ Derive (fun z ⇒ f u z) v) n x y
  end.

Lemma ex_diff_n_ext_loc: ∀ f g n x y,
    locally_2d (fun u v ⇒ f u v = g u v) x y
      → ex_diff_n f n x y → ex_diff_n g n x y.

Lemma ex_diff_n_m :
  ∀ n m, (m ≤ n)%nat → ∀ f x y, ex_diff_n f n x y → ex_diff_n f m x y.

Lemma ex_diff_n_deriv_aux1: ∀ f n x y,
  ex_diff_n f (S n) x y → ex_diff_n (fun u v ⇒ Derive (fun z ⇒ f z v) u) n x y.

Lemma ex_diff_n_deriv_aux2: ∀ f n x y,
  ex_diff_n f (S n) x y → ex_diff_n (fun u v ⇒ Derive (fun z ⇒ f u z) v) n x y.

Lemma ex_diff_n_deriv: ∀ n p q, (p+q ≤ n)%nat → ∀ f x y,
    ex_diff_n f n x y→ ex_diff_n (partial_derive p q f) (n -(p+q)) x y.

Lemma ex_diff_n_ex_deriv_inf_1 : ∀ n p k, (p+k < n)%nat → ∀ f x y,
    ex_diff_n f n x y →
    ex_derive (fun z : R ⇒ partial_derive p k f z y) x.

Lemma ex_diff_n_ex_deriv_inf_2 : ∀ n p k, (p+k < n)%nat → ∀ f x y,
  ex_diff_n f n x y →
  ex_derive (fun z ⇒ partial_derive p k f x z) y.

Lemma ex_diff_n_continuity_inf_1 : ∀ n p k, (p+k < n)%nat → ∀ f x y,
  ex_diff_n f n x y →
  continuity_2d_pt (fun u v ⇒ Derive (fun z : R ⇒ partial_derive p k f z v) u) x y.

Lemma Derive_partial_derive_aux1: ∀ p f x y,
  locally_2d (ex_diff_n f (S p)) x y →
   partial_derive 1 p f x y = partial_derive 0 p (partial_derive 1 0 f) x y.

Lemma Derive_partial_derive_aux2: ∀ p k f x y,
  locally_2d (ex_diff_n f (p+S k)) x y →
  partial_derive 0 1 (partial_derive p k f) x y =
     partial_derive p (S k) f x y.

Lemma Derive_partial_derive: ∀ p k f x y,
  locally_2d (ex_diff_n f (p+S k)) x y →
  Derive (fun v : R ⇒ partial_derive p k f x v) y =
     partial_derive p (S k) f x y.

Lemma ex_diff_n_continuity_inf_2 : ∀ n p k, (p+k < n)%nat → ∀ f x y,
  ex_diff_n f n x y →
  continuity_2d_pt (fun u v ⇒ Derive (fun z : R ⇒ partial_derive p k f u z) v) x y.

Definition DL_regular_n f m x y :=
  ∃ D, locally_2d (fun u v ⇒
    Rabs (f u v - DL_pol m f x y (u-x) (v-y)) ≤ D × (Rmax (Rabs (u-x)) (Rabs (v-y))) ^ (S m)) x y.

Theorem Taylor_Lagrange_2d : ∀ f n x y,
  locally_2d (fun u v ⇒ ex_diff_n f (S n) u v) x y → DL_regular_n f n x y.