forked from coq-community/corn
-
Notifications
You must be signed in to change notification settings - Fork 0
/
CPoly_Lagrange.v
172 lines (152 loc) · 5.19 KB
/
CPoly_Lagrange.v
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
Require Import
Unicode.Utf8
Setoid Arith List Program Permutation
CSetoids CRings CPoly_Degree CPoly_ApZero
CRArith Qmetric Qring CReals
stdlib_omissions.Pair stdlib_omissions.Q
list_separates SetoidPermutation.
Require ne_list.
Import ne_list.notations.
Open Local Scope CR_scope.
Local Notation Σ := cm_Sum.
Local Notation Π := cr_Product.
Section contents.
Notation QPoint := (Q * CR)%type.
Notation CRPoint := (CR * CR)%type.
Variables
(qpoints: list QPoint)
(distinct: QNoDup (map fst qpoints)).
Let crpoints := map (first inject_Q_CR) qpoints.
(** Definition of the Lagrange polynomial: *)
Definition L: cpoly CRasCRing :=
Σ (map (fun p => let '((x, y), rest) := p in
_C_ y [*] Π (map (fun xy' => (' (- fst xy')%Q [+X*] [1]) [*] _C_ (' (/ (x - fst xy')))) rest))
(separates qpoints)).
(** Its degree follows easily from its structure: *)
Lemma degree: degree_le (length (tl qpoints)) L.
Proof with auto.
unfold L.
apply degree_le_Sum. intros ? H.
destruct (proj1 (in_map_iff _ _ _) H) as [[[] v] [[] B]]. clear H.
apply degree_le_mult_constant_l.
clear crpoints.
destruct qpoints.
simpl in B.
exfalso...
simpl length.
replace (@length (prod Q (RegularFunction Q_as_MetricSpace)) l)
with (length (map (fun xi => (' (- fst xi)%Q[+X*][1])[*]_C_ (' (/ (q - fst xi)))) v) * 1)%nat.
apply degree_le_Product.
intros.
destruct (proj1 (in_map_iff _ _ _) H) as [?[[]?]].
apply degree_le_mult_constant_r.
apply degree_le_cpoly_linear_inv.
apply (degree_le_c_ CRasCRing [1]).
rewrite map_length.
ring_simplify.
apply eq_add_S.
rewrite <- (separates_elem_lengths (p0 :: l) v)...
replace v with (snd ((q, s), v))...
apply in_map...
Qed.
(** Applying the polynomial gives what you'd expect: *)
Definition functional (x: Q): CR :=
Σ (map (fun p => let '((xj, y), rest) := p in y [*] ' Π (map (fun xi => (x - fst xi) * / (xj - fst xi))%Q rest))
(separates qpoints)).
Lemma apply x: (L ! ' x) [=] functional x.
Proof.
unfold L, functional.
rewrite cm_Sum_apply.
autorewrite with apply.
rewrite map_map.
apply (@cm_Sum_eq _). intros [[u v] w].
autorewrite with apply.
apply mult_wd. reflexivity.
rewrite inject_Q_product.
rewrite cr_Product_apply.
do 2 rewrite map_map.
apply (@cm_Sum_eq (Build_multCMonoid CRasCRing)).
intros [p q].
rewrite <- CRmult_Qmult.
autorewrite with apply.
change (((' (- p)%Q + ' x * (1 + ' x * 0)) * ' (/ (u - p))) == (' (x + - p)%Q * ' (/ (u - p))))%CR.
rewrite <- CRplus_Qplus.
generalize (/ (u - p)).
intros. ring.
Qed.
(** Finally, the polynomial actually interpolates the given points: *)
Lemma interpolates: interpolates crpoints L.
Proof with auto.
intros xy.
unfold crpoints. rewrite in_map_iff.
intros [[xi y] [V W]].
subst. simpl @fst. simpl @snd.
rewrite apply. unfold functional.
destruct (move_to_front _ _ W) as [x H1].
set (fun p => let '(xj, y1, rest) := p in y1[*] ' Π (map (fun xi0 : Q and CR => ((xi - fst xi0) * / (xj - fst xi0))%Q) rest)).
assert ((pair_rel eq (@Permutation _) ==> @st_eq _) s s) as H2.
intros [u w] [[p q]] [A B].
simpl in A, B. subst u. cbv beta iota.
apply mult_wd. reflexivity.
apply inject_Q_CR_wd.
rewrite B. reflexivity.
pose proof (separates_Proper _ _ H1) as H3.
assert (@Equivalence ((Q * CR) * list (Q * CR)) (pair_rel (@eq _) (@Permutation _))) as T.
apply Pair.Equivalence_instance_0.
apply _.
apply _.
(* I'm really clueless why this rewrite has ever worked? *)
etransitivity. apply cm_Sum_Proper. apply cag_commutes. symmetry.
apply (@map_perm_proper _ CR (pair_rel eq (@Permutation _)) (@st_eq _) T _ _ _ H2 _ _ H3).
clear H2 H3.
subst s.
simpl @cm_Sum.
rewrite cm_Sum_units.
setoid_replace (' Π (map (fun xi0 : Q and CR => ((xi - fst xi0) * / (xi - fst xi0))%Q) x)) with 1%CR.
change (y * 1 + 0 == y). ring.
apply inject_Q_CR_wd.
rewrite cr_Product_ones. reflexivity.
intros.
destruct (proj1 (in_map_iff _ _ _) H) as [x1 [[] H4]].
simpl.
apply Qmult_inv_r.
unfold QNoDup in distinct.
revert distinct.
apply Permutation_sym in H1. (* only needed to work around evar anomaly *)
rewrite H1.
intros H3 H5.
simpl in H3.
inversion_clear H3.
apply H0.
apply in_map_iff.
exists (fst x1).
split...
apply Qred_complete.
symmetry...
apply -> Qminus_eq...
apply in_map...
intros.
destruct (proj1 (in_map_iff _ _ _) H) as [[[v w] u] [[] H4]]. clear H.
destruct (proj1 (in_map_iff _ _ _) H4) as [[[n m] k] [A B]].
inversion_clear A.
simpl @cr_Product.
setoid_replace (xi - xi)%Q with 0%Q by (simpl; ring).
repeat rewrite Qmult_0_l.
change ((w: CR) * 0 == 0).
ring.
Qed. (* Todo: Clean up more. *)
End contents.
Lemma interpolates_economically (qpoints: ne_list (Q * CR)): QNoDup (map fst qpoints) →
interpolates_economically (ne_list.map (first inject_Q_CR) qpoints) (L qpoints).
Proof with auto.
split.
rewrite ne_list.list_map.
apply interpolates...
rewrite ne_list.list_map.
rewrite tl_map.
rewrite map_length.
apply degree...
Qed.
Module notations.
Notation lagrange_poly := L.
End notations.