end of merge

Knuth_axiom5.v

@@ -2,7 +2,6 @@ From mathcomp Require Import ssreflect ssrbool ssrnat eqtype ssrfun seq tuple.
 From mathcomp Require Import choice path finset finfun fintype bigop fingraph.
 From mathcomp Require Import ssralg ssrnum matrix mxalgebra perm fingroup.
 From mathcomp Require Import all_algebra.
-(* From mathcomp Require Import finmap. *)
 From mathcomp Require Import zmodp.
 
 Require Import triangulation_algorithm.
@@ -42,9 +41,9 @@ Proof. by elim/elimI3: i. Qed.
 
 Section abstract_over_types.
 
-Variable P : finType.
+Variable P : Type.
 
-Variable R : realDomainType.
+Variable R : comUnitRingType.
 
 Variable coords : P -> R * R.
 
@@ -52,33 +51,33 @@ Variable p0 : P.
 
 Notation "<< a ; b ; c >>" := (mkf3 a b c).
 
-Notation surface_mx := (surface_mx coords).
-
-Notation ccw := (ccw coords).
-
-Lemma expand_det_surface  a b c :
+Lemma expand_det_surface  (a b c : R * R):
   \det (surface_mx << a; b; c >>) =
-   (coords b).1 * (coords c).2 - (coords b).2 * (coords c).1 -
-    ((coords a).1 * (coords c).2 - (coords a).2 * (coords c).1) +
-    ((coords a).1 * (coords b).2 - (coords a).2 * (coords b).1).
+   b.1 * c.2 - b.2 * c.1 - (a.1 * c.2 - a.2 * c.1) + a.1 * b.2 - a.2 * b.1.
 Proof.
 repeat rewrite (expand_det_col _ ord0) /cofactor /row' /col' !big_ord_recr
    big_ord0 /= add0r !(mxE, ffunE, addn0, expr0, expr1, expr2, mxE, ffunE, 
-   det_mx00, mul1r, mulNr, mulrN, opprK, mulr1) /=.
-by rewrite !(fun f => mulrC (f (coords c)),
-          fun f g => mulrC (f (coords b)) (g (coords a))).
+   det_mx00, mul1r, mulNr, mulrN, opprK, mulr1, addrA) /=.
+by rewrite !(mulrC c.1, mulrC b.1 a.2, mulrC c.2, mulrC b.2 a.1).
 Qed.
 
-Definition R' := (R : Type).
+Definition ptb (t : R) (p : R * R) :=
+  (p.1 + t * p.2, p.2 + t * (p.1 ^ 2 + p.2 ^ 2)).
+
+Definition ptb' '(a, b) : R * R := (a, a ^ 2 + b ^ 2).
 
-Definition mul : R' -> R' -> R' := @GRing.mul _.
-Definition add : R' -> R' -> R' := @GRing.add _.
-Definition sub : R' -> R' -> R' := (fun x y => x - y).
-Definition opp : R' -> R' := @GRing.opp _.
-Definition zero : R' := 0.
-Definition one : R' := 1.
+Definition ptb'' '(a, b) : R * R := (b, a ^ 2 + b ^ 2).
 
-Definition R2_theory :=
+Let R' := (R : Type).
+
+Let mul : R' -> R' -> R' := @GRing.mul _.
+Let add : R' -> R' -> R' := @GRing.add _.
+Let sub : R' -> R' -> R' := (fun x y => x - y).
+Let opp : R' -> R' := @GRing.opp _.
+Let zero : R' := 0.
+Let one : R' := 1.
+
+Let R2_theory :=
   @mk_rt R' zero one add mul sub opp
    (@eq R')
    (@add0r R) (@addrC R) (@addrA R) (@mul1r R) (@mulrC R)
@@ -93,53 +92,176 @@ rewrite ?mxE /= ?(expr0, exprS, mulrS, mulr0n) -?[@GRing.add _]/add
 match goal with |- @eq ?X _ _ => change X with R' end;
 ring.
 
-Lemma Knuth_axiom5 a b c d e :
-  0 < \det (surface_mx << a; b; c >>) ->
-  0 < \det (surface_mx << a; b; d >>) ->
-  0 < \det (surface_mx << a; b; e >>) ->
-  0 < \det (surface_mx << a; c; d >>) ->
-  0 < \det (surface_mx << a; d; e >>) ->
-  0 < \det (surface_mx << a; c; e >>).
+Lemma expand_ptb (t : R) (a b c : R * R):
+\det (surface_mx << ptb t a; ptb t b; ptb t c >>) =
+\det (surface_mx << a; b; c>>) +
+\det (surface_mx << ptb' a; ptb' b; ptb' c >>) * t +
+\det(surface_mx << ptb'' a;  ptb'' b; ptb'' c>>) * t ^ 2.
+Proof.
+rewrite !expand_det_surface /ptb /=.
+rewrite /exprz !exprS !mulr1 /ptb' /ptb''.
+case: a b c=> [a1 a2][b1 b2][c1 c2] /=.
+rewrite /exprz; mc_ring.
+Qed.
+
+(* TODO : This lemma needs renaming. *)
+Lemma step (a b c : R * R) :
+  \det (surface_mx << a; b; c >>) = 0 ->
+  (b.1 - a.1) * c.2 = (b.2 - a.2) * c.1 - a.1 * b.2 + a.2 * b.1.
+Proof.
+rewrite expand_det_surface  opprB !addrA 2!(addrAC _ _ (- (a.1 * c.2))) -mulrBl.
+rewrite -!addrA=> /eqP; rewrite addr_eq0=> /eqP ->.
+mc_ring.
+Qed.
+
+Lemma det_lin3 (a1 a2 b1 b2 c1 c2 v1 v2 v3 v4 : R) :
+  \det (surface_mx << (a1, a2 + v4 * v1); (b1, b2 + v4 * v2);
+          (c1, c2 + v4 * v3) >>) =
+  \det (surface_mx << (a1, a2); (b1, b2); (c1, c2) >>) +
+  v4 * \det (surface_mx << (a1, v1); (b1, v2); (c1, v3) >>).
+Proof.
+rewrite !expand_det_surface /=; mc_ring.
+Qed.
+
+Lemma det_swap23 (a1 a2 b1 b2 c1 c2 : R) :
+ \det (surface_mx << (a1, a2); (b1, b2); (c1, c2) >>) =
+ -\det (surface_mx << (a2, a1); (b2, b1); (c2, c1) >>).
+Proof.
+rewrite !expand_det_surface /=; mc_ring.
+Qed.
+
+Lemma det_D3 (a1 a2 b1 b2 c1 c2 v1 v2 v3 : R) :
+  \det (surface_mx << (a1, a2 + v1); (b1, b2 + v2);
+          (c1, c2 + v3) >>) =
+  \det (surface_mx << (a1, a2); (b1, b2); (c1, c2) >>) +
+  \det (surface_mx << (a1, v1); (b1, v2); (c1, v3) >>).
 Proof.
-suff : \det (surface_mx << a; b; d >>) *
+rewrite !expand_det_surface /=; mc_ring.
+Qed.
+
+Lemma det_M3 (a1 a2 b1 b2 c1 c2 v : R) :
+  \det (surface_mx << (a1, v * a2); (b1, v * b2); (c1, v * c2) >>) =
+  v * \det (surface_mx << (a1, a2); (b1, b2); (c1, c2) >>).
+Proof.
+rewrite !expand_det_surface /=; mc_ring.
+Qed.
+
+Lemma det_013 (a1 b1 c1 : R) :
+  \det (surface_mx << (a1, a1); (b1, b1); (c1, c1) >>) = 0.
+Proof.
+rewrite !expand_det_surface /=; mc_ring.
+Qed.
+
+Lemma det_003 (a1 b1 c1 k : R) :
+  \det (surface_mx << (a1, k); (b1, k); (c1, k) >>) = 0.
+Proof.
+rewrite !expand_det_surface /=; mc_ring.
+Qed.
+
+Lemma expand_van_der_monde (a b c : R) :
+ \det(surface_mx << (a, a ^ 2); (b, b ^ 2); (c, c ^ 2) >>) =
+ (a - b) * (b - c) * (c - a).
+Proof.
+rewrite expand_det_surface /exprz /=; mc_ring.
+Qed.
+
+Lemma expand_ptb' (a b c : R * R) :
+  (b.1 - a.1) * c.2 = (b.2 - a.2) * c.1 - a.1 * b.2 + a.2 * b.1 ->
+  (b.1 - a.1) ^ 2 * \det (surface_mx << ptb' a; ptb' b; ptb' c >>) =
+  ((b.1 - a.1) ^ 2 + (b.2 - a.2) ^ 2) * 
+  (b.1 - a.1) * (c.1 - b.1) * (c.1 - a.1).
+Proof.
+rewrite /ptb'.
+case: a b c=>[a1 a2][b1 b2][c1 c2] /= => cnd.
+rewrite det_D3 mulrDr -[X in _ + X]det_M3.
+have tmp : forall x, (b1 - a1) ^ 2 * x ^ 2 = ((b1 - a1) * x) ^ 2.
+  move=> x; rewrite /exprz; mc_ring.
+rewrite !tmp {tmp}; rewrite cnd.
+have -> : (b1 - a1) * a2 = (b2 - a2) * a1 - a1 * b2 + a2 * b1 by mc_ring.
+have -> : (b1 - a1) * b2 = (b2 - a2) * b1 - a1 * b2 + a2 * b1 by mc_ring.
+have -> : ((b2 - a2) * a1 - a1 * b2 + a2 * b1) ^ 2 =
+          (b2 - a2) ^ 2 * a1 ^ 2 + (2%:R * (b2 - a2) * (a2 * b1 - a1 * b2)) * a1 +
+          (a2 * b1 - a1 * b2) ^ 2.
+  by rewrite /exprz mulr2n; mc_ring.
+have -> : ((b2 - a2) * b1 - a1 * b2 + a2 * b1) ^ 2 =
+          (b2 - a2) ^ 2 * b1 ^ 2 + (2%:R * (b2 - a2) * (a2 * b1 - a1 * b2)) * b1 +
+          (a2 * b1 - a1 * b2) ^ 2.
+  by rewrite /exprz mulr2n; mc_ring.
+have -> : ((b2 - a2) * c1 - a1 * b2 + a2 * b1) ^ 2 =
+          (b2 - a2) ^ 2 * c1 ^ 2 + (2%:R * (b2 - a2) * (a2 * b1 - a1 * b2)) * c1 +
+          (a2 * b1 - a1 * b2) ^ 2.
+  by rewrite /exprz mulr2n; mc_ring.
+rewrite det_D3 det_lin3 det_003 det_M3 det_013 !mulr0 !addr0.
+rewrite -mulrDl -!mulrA; congr (_ * _).
+rewrite expand_van_der_monde; mc_ring.
+Qed.
+
+Lemma expand_ptb'' (a b c : R * R) :
+  (b.1 - a.1) * c.2 = (b.2 - a.2) * c.1 - a.1 * b.2 + a.2 * b.1 ->
+  (b.2 - a.2) ^ 2 * \det (surface_mx << ptb'' a; ptb'' b; ptb'' c >>) =
+  ((b.1 - a.1) ^ 2 + (b.2 - a.2) ^ 2) * 
+  (b.2 - a.2) * (c.2 - b.2) * (c.2 - a.2).
+Proof.
+rewrite /ptb''.
+case: a b c=>[a1 a2][b1 b2][c1 c2] /= => cnd.
+have cnd' : (b2 - a2) * c1 = (b1 - a1) * c2 + a1 * b2 - a2 * b1.
+  by rewrite cnd; mc_ring.
+rewrite det_D3 mulrDr -[X in X + _]det_M3.
+have tmp : forall x, (b2 - a2) ^ 2 * x ^ 2 = ((b2 - a2) * x) ^ 2.
+  move=> x; rewrite /exprz; mc_ring.
+rewrite !tmp {tmp}; rewrite cnd'.
+have -> : (b2 - a2) * a1 = (b1 - a1) * a2 + a1 * b2 - a2 * b1 by mc_ring.
+have -> : (b2 - a2) * b1 = (b1 - a1) * b2 + a1 * b2 - a2 * b1 by mc_ring.
+have -> : ((b1 - a1) * a2 + a1 * b2 - a2 * b1) ^ 2 =
+          (b1 - a1) ^ 2 * a2 ^ 2 + (2%:R * (a1 - b1) * (a2 * b1 - a1 * b2)) * a2 +
+          (a2 * b1 - a1 * b2) ^ 2.
+  by rewrite /exprz mulr2n; mc_ring.
+have -> : ((b1 - a1) * b2 + a1 * b2 - a2 * b1) ^ 2 =
+          (b1 - a1) ^ 2 * b2 ^ 2 + (2%:R * (a1 - b1) * (a2 * b1 - a1 * b2)) * b2 +
+          (a2 * b1 - a1 * b2) ^ 2.
+  by rewrite /exprz mulr2n; mc_ring.
+have -> : ((b1 - a1) * c2 + a1 * b2 - a2 * b1) ^ 2 =
+          (b1 - a1) ^ 2 * c2 ^ 2 + (2%:R * (a1 - b1) * (a2 * b1 - a1 * b2)) * c2 +
+          (a2 * b1 - a1 * b2) ^ 2.
+  by rewrite /exprz mulr2n; mc_ring.
+rewrite det_D3 det_lin3 det_003 det_M3 det_013 !mulr0 !addr0.
+rewrite -mulrDl -!mulrA; congr (_ * _).
+rewrite expand_van_der_monde; mc_ring.
+Qed.
+
+Lemma cramer (a b c d e : R * R) :
+       \det (surface_mx << a; b; d >>) *
        \det (surface_mx << a; c; e >>) =
        \det (surface_mx << a; b; c >>) *
        \det (surface_mx << a; d; e >>) +
        \det (surface_mx << a; b; e >>) *
        \det (surface_mx << a; c; d >>).
-  move => cramer abc abd abe acd ade.
-  by rewrite -(pmulr_lgt0 _ abd) mulrC cramer addr_gt0 // mulr_gt0.
+Proof.
 rewrite !expand_det_surface.
 mc_ring.
 Qed.
 
-Lemma Knuth_axiom5c a b c d e :
-  ccw a b c -> ccw a b d -> ccw a b e -> ccw a c d -> ccw a d e -> ccw a c e.
-Proof. exact: Knuth_axiom5. Qed.
-
-Lemma Knuth_axiom5' a b c d e :
-  0 < \det (surface_mx << a; b; c >>) ->
-  0 < \det (surface_mx << a; b; d >>) ->
-  0 < \det (surface_mx << a; b; e >>) ->
-  0 < \det (surface_mx << b; c; d >>) ->
-  0 < \det (surface_mx << b; d; e >>) ->
-  0 < \det (surface_mx << b; c; e >>).
-Proof.
-suff : \det (surface_mx << a; b; d >>) *
+Lemma cramer' (a b c d e : R * R) :
+       \det (surface_mx << a; b; d >>) *
        \det (surface_mx << b; c; e >>) =
        \det (surface_mx << a; b; c >>) *
        \det (surface_mx << b; d; e >>) +
        \det (surface_mx << a; b; e >>) *
        \det (surface_mx << b; c; d >>).
-  move => cramer abc abd abe acd ade.
-  by rewrite -(pmulr_lgt0 _ abd) mulrC cramer addr_gt0 // mulr_gt0.
+Proof.
 rewrite !expand_det_surface.
 mc_ring.
 Qed.
 
-Lemma Knuth_axiom5c' a b c d e :
-  ccw a b c -> ccw a b d -> ccw a b e -> ccw b c d -> ccw b d e -> ccw b c e.
-Proof. exact: Knuth_axiom5'. Qed.
+Lemma Knuth_4_main (a b c d : R * R) :
+  \det (surface_mx << b; c; d >>) -
+  \det (surface_mx << a; c; d >>) +
+  \det (surface_mx << a; b; d >>) -
+  \det (surface_mx << a; b; c >>) = 0.
+Proof.
+rewrite !expand_det_surface.
+mc_ring.
+Qed.
 
 Definition shift_cycle_val := [ffun i : 'I_3 => 1 + i].
 
@@ -148,7 +270,7 @@ Proof using. by apply/perm_proof/addrI. Qed.
 
 Open Scope group_scope.
 
-Lemma cycle_det_surface1 a b c :
+Lemma cycle_det_surface1 (a b c : R * R) :
   \det (surface_mx << b; c; a >>) = \det (surface_mx << a; b; c>>).
 Proof using.
 have shift_cycle_even : (Perm shift_cycle_proof : bool) = false.
@@ -164,7 +286,7 @@ rewrite -row_permE; apply/matrixP => k l; rewrite !mxE !ffunE.
 by elim/elimI3: k => //=; rewrite unlock ffunE.
 Qed.
 
-Lemma cycle_det_surface f (i : 'I_3) :
+Lemma cycle_det_surface (f : (R * R)%type ^ 3) (i : 'I_3) :
   \det (surface_mx <<f i; f (i + 1%R); f (i - 1%R)>>) =
   \det (surface_mx <<f 0; f 1%R; f (-1%R)>>).
 Proof using.
@@ -175,7 +297,7 @@ Qed.
 
 Close Scope group_scope.
 
-Lemma swap_det_surface a b c :
+Lemma swap_det_surface (a b c : R * R) :
   \det (surface_mx << a; b; c >>) = - \det (surface_mx << b; a; c >>).
 Proof using.
 have m_eq U V W : surface_mx << V; U; W >> =
@@ -185,40 +307,75 @@ have m_eq U V W : surface_mx << V; U; W >> =
 by rewrite m_eq detM det_perm odd_tperm /= expr1 mulNr mul1r.
 Qed.
 
+End abstract_over_types.
+
+Section abstract_over_types'.
+
+Variable P : Type.
+
+Variable R : realDomainType.
+
+Variable coords : P -> R * R.
+
+Variable p0 : P.
+
+Notation "<< a ; b ; c >>" := (mkf3 a b c).
+
+Lemma Knuth_axiom5 (a b c d e : R * R) :
+  0 < \det (surface_mx << a; b; c >>) ->
+  0 < \det (surface_mx << a; b; d >>) ->
+  0 < \det (surface_mx << a; b; e >>) ->
+  0 < \det (surface_mx << a; c; d >>) ->
+  0 < \det (surface_mx << a; d; e >>) ->
+  0 < \det (surface_mx << a; c; e >>).
+Proof.
+move => abc abd abe acd ade.
+by rewrite -(pmulr_lgt0 _ abd) mulrC cramer addr_gt0 // mulr_gt0.
+Qed.
+
+Lemma Knuth_axiom5c a b c d e :
+  ccw coords a b c -> ccw coords a b d -> ccw coords a b e ->
+  ccw coords a c d -> ccw coords a d e -> ccw coords a c e.
+Proof. exact: Knuth_axiom5. Qed.
+
+Lemma Knuth_axiom5' (a b c d e : R * R) :
+  0 < \det (surface_mx << a; b; c >>) ->
+  0 < \det (surface_mx << a; b; d >>) ->
+  0 < \det (surface_mx << a; b; e >>) ->
+  0 < \det (surface_mx << b; c; d >>) ->
+  0 < \det (surface_mx << b; d; e >>) ->
+  0 < \det (surface_mx << b; c; e >>).
+Proof.
+move => abc abd abe acd ade.
+by rewrite -(pmulr_lgt0 _ abd) mulrC cramer' addr_gt0 // mulr_gt0.
+Qed.
+
+Lemma Knuth_axiom5c' a b c d e :
+  ccw coords a b c -> ccw coords a b d -> ccw coords a b e ->
+  ccw coords b c d -> ccw coords b d e -> ccw coords b c e.
+Proof. exact: Knuth_axiom5'. Qed.
+
 Definition fina := (addrN, addNr, add0r).
 
 Lemma Knuth_4 d a b c:
-  ccw a b d -> ccw b c d -> ccw c a d -> ccw a b c.
+  ccw coords a b d -> ccw coords b c d -> ccw coords c a d -> ccw coords a b c.
 Proof using.
-rewrite /ccw.
-suff main : \det(surface_mx << b; c; d >>) -
-       \det(surface_mx << a; c; d >>) +
-       \det(surface_mx << a; b; d >>) -
-       \det(surface_mx << a; b; c >>) = 0.
-  move=> c1 c2 c3; rewrite -oppr0 -main -[X in _ < X]add0r !opprD !opprK.
-  rewrite ltr_add2r subr_lt0 (ltr_trans _ c1) // addrC subr_lt0.
-  by rewrite (ltr_trans _ c2) //; rewrite swap_det_surface oppr_lt0.
-rewrite !expand_det_surface !opprD !addrA !opprK.
-set a1 := (coords a).1; set a2 := (coords a).2.
-set b1 := (coords b).1; set b2 := (coords b).2.
-set c1 := (coords c).1; set c2 := (coords c).2.
-set d1 := (coords d).1; set d2 := (coords d).2.
-rewrite -!(addrAC _ (- (c1 * d2)))  !fina -!(addrAC _ (c2 * d1)) !fina.
-rewrite -!(addrAC _ (b1 * d2)) !fina -!(addrAC _ (- (b2 * d1))) !fina.
-rewrite -!(addrAC _ (- (b1 * c2))) !fina -!(addrAC _ (b2 * c1)) !fina.
-rewrite -!(addrAC _ (- (a1 * d2))) !fina -!(addrAC _ (a2 * d1)) !fina.
-rewrite -!(addrAC _ (a1 * c2)) !fina -!(addrAC _ (- (a2 * c1))) !fina.
-by rewrite -!(addrAC _ (- (a1 * b2))) !fina.
+rewrite /ccw /ccw0.
+have main:= Knuth_4_main (coords a) (coords b) (coords c) (coords d).
+move=> c1 c2 c3; rewrite -oppr0 -main -[X in _ < X]add0r !opprD !opprK.
+rewrite ltr_add2r subr_lt0 (ltr_trans _ c1) // addrC subr_lt0.
+by rewrite (ltr_trans _ c2) //; rewrite swap_det_surface oppr_lt0.
 Qed.
 
-Lemma Knuth_1 a b c : ccw a b c = ccw c a b.
+Lemma Knuth_1 a b c : ccw coords a b c = ccw coords c a b.
 Proof.
-by rewrite /ccw cycle_det_surface1.
+by rewrite /ccw /ccw0 cycle_det_surface1.
 Qed.
 
-Lemma Knuth_2 a b c : ccw a b c -> ~~ccw b a c.
+Lemma Knuth_2 a b c : ccw coords a b c -> ~~ccw coords b a c.
 Proof.
-by rewrite /ccw swap_det_surface oppr_gt0 ltrNge ler_eqVlt negb_or => /andP[].
+rewrite /ccw /ccw0 swap_det_surface oppr_gt0 ltrNge ler_eqVlt negb_or.
+by move => /andP[].
 Qed.
 
-End abstract_over_types.
+End abstract_over_types'.
\ No newline at end of file

runner.v

+From mathcomp Require Import all_ssreflect ssralg ssrnum ssrint.
+
+From triangles Require Import triangles3 to_ccw_system.
+
+Open Scope ring_scope.
+
+Section fix_fin_type.
+
+Variable R : realDomainType.
+
+Variable l : seq (R * R).
+
+Let P := seq_sub l.
+
+Definition no_det_ccf : cc_fun [finType of P] :=
+  @Build_cc_fun _
+  (fun a' b' c' : P =>
+    let a := val a' in let b := val b' in let c := val c' in
+    let dt :=
+        b.1 * c.2 - b.2 * c.1 - (a.1 * c.2 - a.2 * c.1) + a.1 * b.2
+           - a.2 * b.1 in
+     if 0 < dt then true
+     else if 0 == dt then
+       if a.1 < b.1 then (b.1 < c.1) || (c.1 < a.1)
+       else if b.1 < a.1 then (b.1 < c.1) && (c.1 < a.1)
+       else if a.2 < b.2 then (b.2 < c.2) || (c.2 < a.2)
+       else if b.2 < a.2 then (b.2 < c.2) && (c.2 < a.2)
+       else false
+    else false).
+
+Lemma no_det_ccfP (a b c : P) :
+  triangles3.ccw no_det_ccf a b c = ccw (val a) (val b) (val c).
+Proof.
+by rewrite /ccw /= /surface Knuth_axiom5.expand_det_surface.
+Qed.
+
+End fix_fin_type.
+
+Section assoc_map.
+
+Variable P : eqType.
+
+Variable default : P.
+
+Definition lsort := list (P * P).
+
+Definition lmap0 : lsort := nil.
+
+Definition map_size : lsort -> nat := 
+  fun m => (2 * (size m) + 3)%N.
+
+Definition update : lsort -> P -> P -> lsort :=
+  fun l k v =>
+    (fix rec l :=
+      match l with 
+      | nil => (k, v) :: nil
+      | (k', v')::tl => if k == k' then (k', v)::tl else (k', v')::rec tl
+      end) l.
+
+Definition eval : lsort -> P -> P :=
+  fun l k =>
+    (fix rec l :=
+        match l with
+        | nil => default
+        | (k', v')::tl => if k == k' then v' else rec tl
+        end) l.
+
+Lemma leval_update1 : forall (m : lsort) x y, eval (update m x y) x = y.
+Proof.
+elim=> [ | [k v] tl IH] x y /=; case: ifP=> //.
+    by rewrite eqxx.
+  by move=> /eqP -> /=; rewrite eqxx.
+by rewrite /= => ->.
+Qed.
+
+Lemma leval_update_dif : forall (m : lsort) x x' y, x != x' ->
+    eval (update m x y) x' = eval m x'.
+Proof.
+elim=> [ | [k v] tl IH] x x' y /=; case: ifP=> //.
+    by rewrite eq_sym => ->.
+  by move=> /eqP -> /=; rewrite eq_sym => /negbTE ->.
+by rewrite /= => xnk xn'; case:ifP=> _ //; apply: IH.
+Qed.
+
+Lemma map_n_default_pair (m : lsort) (x : P) :
+  eval m x != default -> (x, (eval m x)) \in m.
+Proof.
+elim: m=> [|[x' v'] tl IH]; first by rewrite eqxx.
+rewrite /=; case: ifP; rewrite !inE; first by move=> /eqP ->; rewrite eqxx.
+by move=> xnx' vndef; rewrite IH ?orbT.
+Qed.
+
+Lemma lmap_sizeP (m : lsort) (x : P) (p : seq P):
+                (fpath (eval m)) x p ->
+                uniq (x :: p) ->
+                (size (x :: p) <= map_size m)%N .
+Proof.
+move: p x.
+have main p x : default \notin (x :: p) -> fpath (eval m) x p ->
+   uniq (x :: p) -> (size (x :: p) <= (size m).+1)%N.
+  move=> dnp pp up /=; rewrite ltnS.
+  suff h : {subset p <= [seq i.2 | i <- m]}.
+    have -> : size m = size [seq i.2 | i <- m] by rewrite size_map.
+    by apply: uniq_leq_size => //; move: up; rewrite cons_uniq=> /andP[].
+  move/pathP: pp => /= pp z /(nthP x) [i /= ilts <-] {z}.
+  rewrite -(eqP (pp x i ilts)).
+  have /map_n_default_pair it : eval m (nth x (x :: p) i) != default.
+    rewrite (eqP (pp _ _ ilts)); apply/eqP=> abs; case/negP: dnp; rewrite -abs.
+    by rewrite inE mem_nth ?orbT.
+  by apply/mapP; exists (nth x (x :: p) i, eval m (nth x (x :: p) i)).
+move=> p x; have [] := boolP (default \in x :: p); last first.
+  move=> dnin pp up; move: (main _ _ dnin pp up) => pl.
+  apply: (leq_trans pl); rewrite /map_size addnS ltnS mulSn mul1n.
+  by rewrite -addnA leq_addr.
+move=> din pp up; have/(nth_index x):= din; set i := index _ _.
+have /= : (i < size (x :: p))%N by rewrite /i index_mem.
+rewrite ltnS.
+case: i => [ | i] /= i_s.
+  move=> xdef; move: pp up; rewrite xdef.
+  case: p din i_s => [ | y p] din _; first by rewrite /map_size addnS.
+  move=> /= /andP [ed pp] /andP [dn up].
+  have h : (size (y :: p) <= (size m).+1)%N by apply: main.
+  rewrite /= /map_size mulSn mul1n addnS ltnS (leq_trans h) //.
+  by rewrite addnS ltnS -!addnA leq_addr.
+move=> idef.
+have := cat_take_drop i.+1 (x :: p).
+set p1 := take _ _; set p2 := drop _ _ => xpq.
+have p2shape : p2 = default :: drop i.+1 p by rewrite -idef -drop_nth.
+move: (xpq); rewrite /p1 /= => [[]] pq.
+move: up; rewrite -xpq cat_uniq=> /andP[up1 /andP[nhas up2]].
+have sp1 : (size p1 <= (size m).+1)%N.
+rewrite /p1 -[take _ _]/(x :: take i p); apply: main=> //.
+    by move: nhas; rewrite p2shape /= negb_or /p1 /==> /andP [].
+  by move: (pp); rewrite -{1}pq cat_path => /andP[].
+have sp2 : (size p2 <= (size m).+2)%N.
+  rewrite p2shape; case cd : (drop i.+1 p) => [// | z p'].
+  rewrite -[size [:: _, _ & _]]/(size (z :: p')).+1 ltnS.
+  have uzp : uniq (z :: p').
+    by move: up2; rewrite p2shape cons_uniq cd=> /andP[].
+  apply: main => //.
+    by move: up2; rewrite p2shape cons_uniq cd => /andP[].
+  move: pp; rewrite -pq cat_path p2shape cd /=.
+  by move=> /andP[] ? /andP[] ? /andP[] ?.
+rewrite -pq size_cat -addSn -[S _]/(size p1) /map_size mulSn mul1n.
+rewrite !addnS addn0 -addSn -2!addnS.
+by apply: leq_add.    
+Qed.
+
+Definition assoc_list_map_funs : map_funs P :=
+  @Build_map_funs P lsort nil eval update map_size.
+
+Definition assoc_list_map_system : map_system assoc_list_map_funs.
+constructor.
+  exact leval_update1.
+  exact leval_update_dif.
+  exact lmap_sizeP.
+Qed.
+
+End assoc_map.
+
+Lemma inl0 (T : eqType) (x : T) (l : seq T) :
+  x \in x :: l.
+Proof.
+by rewrite inE eqxx.
+Qed.
+
+Definition map_to_orbit (T : eqType) (x : T) (mf : map_funs T) (m : sort mf):
+  seq T :=
+  x ::(fix aux (n : nat) (y stop : T) :=
+   match n with
+     0 => [::]
+   | S p => if y == stop then [::] else y::aux p (triangles3.eval m y) stop end)
+     (triangles3.map_size m) (triangles3.eval m x) x.
+
+Definition int_naive5 :
+ seq (int * int) -> (seq (seq (int * int)) * seq (int * int) * (int * int)):=
+ fun l =>
+   let l' := (0:int, 0:int) :: l in
+   let p0 := (SeqSub (inl0 _ (0:int, 0:int) l)) in
+   let lp := [seq (insubd p0 x) | x <- l] : seq (seq_sub l') in
+   let ccf := no_det_ccf _ l' in
+   let mf := assoc_list_map_funs [eqType of (seq_sub l')] p0 in
+   let result :=
+     naive5 p0 ccf mf lp in
+     ([seq [seq val x | x <- t] | t <- result.1.1], 
+      [seq val v | v <- map_to_orbit _ result.2 mf result.1.2],
+      (val result.2)).
+

triangulation_algorithm.v

@@ -8,11 +8,46 @@ Set Implicit Arguments.
 Unset Strict Implicit.
 Unset Printing Implicit Defensive.
 
+Section num_computations.
+
+Variable P : Type.
+
+Variable R : unitRingType.
+
+Variables coords : P -> R * R.
+
+Definition mkf3 (T : Type) (a b c : T) :=
+  [ffun i : 'I_3 => if val i == 0 then a else if val i == 1 then b else c].
+
+Notation "<< a ; b ; c >>" := (mkf3 a b c).
+
+Definition project_p (p : R * R) : R ^ 3 :=
+  [ffun i => if val i == 0 then 1%R else if val i == 1 then p.1 else p.2].
+
+Definition surface_mx (t : (R * R)%type ^ 3) :=
+  (\matrix_(i < 3, j < 3) project_p (t i) j)%R.
+
+End num_computations.
+
+Section num_computations_comparison.
+
+Variable P : Type.
+
+Variable R : realDomainType.
+
+Variable coords : P -> R * R.
+
+Definition ccw0 (a b c : R * R) := (0 < \det (surface_mx (mkf3 a b c)))%R.
+
+Definition ccw (a b c : P) := ccw0 (coords a) (coords b) (coords c).
+
+End num_computations_comparison.
+
 Section abstract_over_types.
 
 Variable P : finType.
 
-Variable R : numDomainType.
+Variable R : realDomainType.
 
 Variable coords : P -> R * R.
 
@@ -22,20 +57,9 @@ Variable pick_set : {set P} -> option P.
 
 Variable pick_triangle : {set {set P}} -> P -> option {set P}.
 
-Definition mkf3 (T : Type) (a b c : T) :=
-  [ffun i : 'I_3 => if val i == 0 then a else if val i == 1 then b else c].
-
 Notation "<< a ; b ; c >>" := (mkf3 a b c).
 
-Definition project_p (p : P) : R ^ 3 :=
-  [ffun i =>
-     if val i == 0 then 1%R
-     else if val i == 1 then (coords p).1 else (coords p).2].
-
-Definition surface_mx (t : P ^ 3) :=
-  (\matrix_(i < 3, j < 3) project_p (t i) j)%R.
-
-Definition ccw (a b c : P) := (0 < \det (surface_mx << a; b; c>>))%R.
+Notation ccw := (ccw coords).
 
 Definition separated (t : {set P}) (q p : P) :=
   (q \in t) &&
@@ -280,15 +304,15 @@ Definition naive_triangulate' (s : {set P}) := naive_triangulate #|s| s.
 Definition flip_edge (t1 t2 : {set P}) :=
   [set x |: (t1 :\: t2) :|: (t2 :\: t1) | x in t1 :&: t2].
 
-Definition in_circle_l (p : P) : R ^ 4 :=
+Definition in_circle_l (p : R * R) : R ^ 4 :=
   [ffun i : 'I_4 =>
    if val i < 3 then
      project_p p (inord i)
    else
-    ((coords p).1 ^+ 2 + (coords p).2 ^+ 2)%R].
+    (p.1 ^+ 2 + p.2 ^+ 2)%R].
 
 Definition in_circle_mx (t : P ^ 4) :=
-  (\matrix_(i < 4, j < 4) in_circle_l (t i) j)%R.
+  (\matrix_(i < 4, j < 4) in_circle_l (coords (t i)) j)%R.