diff --git a/Makefile b/Makefile
index b3f0d6162a..3e65f9f3c7 100644
--- a/Makefile
+++ b/Makefile
@@ -88,7 +88,9 @@ EXCLUDED_VOFILES := $(filter $(EXCLUDE_PATTERN),$(VOFILES))
 # add files to this list to prevent them from being built as final
 # targets by the "lite" target
 LITE_UNMADE_VOFILES := src/Curves/Weierstrass/AffineProofs.vo \
-	src/Curves/Weierstrass/Jacobian.vo \
+	src/Curves/Weierstrass/Jacobian/Jacobian.vo \
+	src/Curves/Weierstrass/Jacobian/CoZ.vo \
+	src/Curves/Weierstrass/Jacobian/ScalarMult.vo \
 	src/Curves/Weierstrass/Projective.vo \
 	src/Rewriter/RulesGood.vo \
 	src/Rewriter/All.vo \
@@ -96,7 +98,9 @@ LITE_UNMADE_VOFILES := src/Curves/Weierstrass/AffineProofs.vo \
 	$(EXCLUDED_VO)
 NOBIGMEM_UNMADE_VOFILES := \
 	src/Curves/Weierstrass/AffineProofs.vo \
-	src/Curves/Weierstrass/Jacobian.vo \
+	src/Curves/Weierstrass/Jacobian/Jacobian.vo \
+	src/Curves/Weierstrass/Jacobian/CoZ.vo \
+	src/Curves/Weierstrass/Jacobian/ScalarMult.vo \
 	src/Curves/Weierstrass/Projective.vo \
 	$(PERFTESTING_VO) \
 	$(EXCLUDED_VO)
diff --git a/src/Curves/Weierstrass/Jacobian/CoZ.v b/src/Curves/Weierstrass/Jacobian/CoZ.v
new file mode 100644
index 0000000000..9925cade73
--- /dev/null
+++ b/src/Curves/Weierstrass/Jacobian/CoZ.v
@@ -0,0 +1,872 @@
+Require Import Coq.Classes.Morphisms.
+
+Require Import Crypto.Spec.WeierstrassCurve.
+Require Import Crypto.Curves.Weierstrass.Affine Crypto.Curves.Weierstrass.AffineProofs.
+Require Import Crypto.Curves.Weierstrass.Jacobian.Jacobian.
+Require Import Crypto.Util.Decidable Crypto.Algebra.Field.
+Require Import Crypto.Util.Tactics.BreakMatch.
+Require Import Crypto.Util.Tactics.DestructHead.
+Require Import Crypto.Util.Tactics.SpecializeBy.
+Require Import Crypto.Util.Tactics.SetoidSubst.
+Require Import Crypto.Util.Notations Crypto.Util.LetIn.
+Require Import Crypto.Util.Sum Crypto.Util.Prod Crypto.Util.Sigma.
+Require Import Crypto.Util.FsatzAutoLemmas.
+Require Import Coq.PArith.BinPos.
+
+Module Jacobian.
+  Section Co_Z.
+    Context {F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv} {a b:F}
+            {field:@Algebra.Hierarchy.field F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv}
+            {char_ge_3:@Ring.char_ge F Feq Fzero Fone Fopp Fadd Fsub Fmul (BinNat.N.succ_pos (BinNat.N.two))}
+            {Feq_dec:DecidableRel Feq}.
+    Local Infix "=" := Feq : type_scope. Local Notation "a <> b" := (not (a = b)) : type_scope.
+    Local Notation "0" := Fzero. Local Notation "1" := Fone.
+    Local Infix "+" := Fadd. Local Infix "-" := Fsub.
+    Local Infix "*" := Fmul. Local Infix "/" := Fdiv.
+    Local Notation "x ^ 2" := (x*x). Local Notation "x ^ 3" := (x^2*x).
+    Local Notation "2" := (1+1). Local Notation "4" := (2+2).
+    Local Notation "8" := (4+4). Local Notation "27" := (4*4 + 4+4 +1+1+1).
+    Local Notation Wpoint := (@W.point F Feq Fadd Fmul a b).
+    Context {char_ge_12:@Ring.char_ge F Feq Fzero Fone Fopp Fadd Fsub Fmul 12%positive}
+            {discriminant_nonzero:id(4*a*a*a + 27*b*b <> 0)}.
+
+    Local Notation Wopp := (@W.opp F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation Wadd := (@W.add F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv field Feq_dec char_ge_3 a b).
+    Local Notation Wzero := (@W.zero F Feq Fadd Fmul a b).
+    Local Notation Weq := (@W.eq F Feq Fadd Fmul a b).
+
+    Local Notation point := (@Jacobian.point F Feq Fzero Fadd Fmul a b Feq_dec).
+    Local Notation eq := (@Jacobian.eq F Feq Fzero Fadd Fmul a b Feq_dec).
+    Local Notation opp := (@Jacobian.opp F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation add := (@Jacobian.add F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field char_ge_3).
+    Local Notation double := (@Jacobian.double F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation of_affine := (@Jacobian.of_affine F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation to_affine := (@Jacobian.to_affine F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+
+   (* Stolen from Jacobian.v *)
+    Ltac prept_step_opt :=
+      match goal with
+      | [ H : ?T |- ?T ] => exact H
+      | [ |- ?x = ?x ] => reflexivity
+      | [ H : ?T, H' : ~?T |- _ ] => solve [ exfalso; apply H', H ]
+      end.
+
+    Ltac prept_step :=
+      match goal with
+      | _ => progress prept_step_opt
+      | _ => progress intros
+      (*| _ => progress specialize_by_assumption*)
+      (*| _ => progress specialize_by trivial*)
+      | _ => progress cbv [proj1_sig fst snd] in *
+      | _ => progress autounfold with points_as_coordinates in *
+      | _ => progress destruct_head'_True
+      | _ => progress destruct_head'_unit
+      | _ => progress destruct_head'_prod
+      | _ => progress destruct_head'_sig
+      | _ => progress destruct_head'_and
+      | _ => progress destruct_head'_sum
+      | _ => progress destruct_head'_bool
+      | _ => progress destruct_head'_or
+      | H: context[@dec ?P ?pf] |- _ => destruct (@dec P pf)
+      | |- context[@dec ?P ?pf]      => destruct (@dec P pf)
+      | |- ?P => lazymatch type of P with Prop => split end
+      end.
+    Ltac prept := repeat prept_step.
+    Ltac t := prept; trivial; try contradiction; fsatz.
+
+    Definition x_of (P : point) : F :=
+      match proj1_sig P with
+      | (x, y, z) => x
+      end.
+
+    Definition y_of (P : point) : F :=
+      match proj1_sig P with
+      | (x, y, z) => y
+      end.
+
+    Definition z_of (P : point) : F :=
+      match proj1_sig P with
+      | (x, y, z) => z
+      end.
+
+    Create HintDb points_as_coordinates discriminated.
+    Hint Unfold Proper respectful Reflexive Symmetric Transitive : points_as_coordinates.
+    Hint Unfold Jacobian.point Jacobian.eq W.eq W.point W.coordinates proj1_sig fst snd : points_as_coordinates.
+    Hint Unfold x_of y_of z_of : points_as_coordinates.
+
+        Local Ltac clear_neq :=
+      repeat match goal with
+             | [ H : _ <> _ |- _ ] => clear H
+             end.
+    Local Ltac clear_eq :=
+      repeat match goal with
+             | [ H : _ = _ |- _ ] => clear H
+             end.
+    Local Ltac clear_eq_and_neq :=
+      repeat (clear_eq || clear_neq).
+    Local Ltac clean_up_speed_up_fsatz :=
+      repeat match goal with
+             | [ H : forall a : F, _ = _ -> _ = _ |- _ ] => clear H
+             | [ H : forall a : F, _ <> _ -> _ = _ |- _ ] => clear H
+             | [ H : forall a : F, _ <> _ -> _ <> _ |- _ ] => clear H
+             | [ H : forall a b : F, _ = _ |- _ ] => clear H
+             | [ H : forall a b : F, _ = _ -> _ = _ |- _ ] => clear H
+             | [ H : forall a b : F, _ <> _ -> _ = _ |- _ ] => clear H
+             | [ H : forall a b : F, _ <> _ -> _ <> _ |- _ ] => clear H
+             | [ H : forall a b : F, _ = _ -> _ <> _ -> _ = _ |- _ ] => clear H
+             | [ H : forall a b : F, _ <> _ -> _ <> _ -> _ = _ |- _ ] => clear H
+             | [ H : forall a b : F, _ <> _ -> _ <> _ -> _ <> _ |- _ ] => clear H
+             | [ H : forall a b c : F, _ |- _ ] => clear H
+             | [ H : ?T, H' : ?T |- _ ] => clear H'
+             | [ H : ?x = ?x |- _ ] => clear H
+             | [ H : ?x <> 0, H' : ?x = 1 |- _ ] => clear H
+             | [ H : ?x = 0 |- _ ] => is_var x; clear H x
+             | [ H : ?x = 1 |- _ ] => is_var x; clear H x
+             | [ H : Finv (?x^3) <> 0, H' : ?x <> 0 |- _ ] => clear H
+             end.
+    Local Ltac find_and_do_or_prove_by ty by_tac tac :=
+      lazymatch goal with
+      | [ H' : ty |- _ ]
+        => tac H'
+      | _
+        => assert ty by by_tac
+      end.
+    Local Ltac find_and_do_or_prove_by_fsatz ty tac :=
+      find_and_do_or_prove_by ty ltac:(clear_eq_and_neq; intros; fsatz) tac.
+    Local Ltac find_and_do_or_prove_by_fsatz_with_neq ty tac :=
+      find_and_do_or_prove_by ty ltac:(clear_neq; intros; fsatz) tac.
+    Local Ltac find_and_goal_apply_or_prove_by_fsatz ty :=
+      find_and_do_or_prove_by_fsatz ty ltac:(fun H' => apply H').
+    Local Ltac find_and_apply_or_prove_by_fsatz H ty :=
+      find_and_do_or_prove_by_fsatz ty ltac:(fun H' => apply H' in H).
+    Local Ltac find_and_apply_or_prove_by_fsatz2 H0 H1 ty :=
+      find_and_do_or_prove_by_fsatz ty ltac:(fun H' => apply H' in H0; [ | exact H1 ]).
+    Local Ltac find_and_apply_or_prove_by_fsatz' H ty preapp :=
+      find_and_do_or_prove_by_fsatz ty ltac:(fun H' => let H' := preapp H' in apply H' in H).
+    Local Ltac speed_up_fsatz_step_on pkg :=
+      let goal_if_var_neq_zero (* we want to keep lazymatch below, so we hack backtracking on is_var *)
+          := match goal with
+             | [ |- ?x <> 0 ] => let test := match goal with _ => is_var x end in
+                                 constr:(x <> 0)
+             | _ => I
+             end in
+      lazymatch goal with
+      | [ H : False |- _ ] => solve [ exfalso; exact H ]
+      | [ H : ?T |- ?T ] => exact H
+      | [ H : ?x <> ?x |- _ ] => solve [ exfalso; apply H; reflexivity ]
+      | [ H : ?x = ?y, H' : ?x <> ?y |- _ ] => solve [ exfalso; apply H', H ]
+      | [ H : ?x = 0 |- ?x * ?y <> 0 ] => exfalso
+      | [ H : ?y = 0 |- ?x * ?y <> 0 ] => exfalso
+      | [ H : ~?T |- ?T ] => exfalso
+      | [ H : ?T |- ~?T ] => exfalso
+      | [ H : ?x - ?x = 0 |- _ ] => clear H
+      | [ H : ?x * ?y = 0, H' : ?x = 0 |- _ ] => clear H
+      | [ H : ?x * ?y = 0, H' : ?y = 0 |- _ ] => clear H
+      | [ |- goal_if_var_neq_zero ] => intro
+      | [ H : ?x = Fopp ?x |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a, a = Fopp a -> a = 0)
+      | [ H : ?x <> Fopp ?x |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a, a <> Fopp a -> a <> 0)
+      | [ H : ?x + ?x = 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall y, y + y = 0 -> y = 0)
+      | [ H : Finv (?x^3) = 0, H' : ?x <> 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a, Finv (a^3) = 0 -> a <> 0 -> False)
+      | [ H : ?x - ?y = 0, H' : ?y <> ?x |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a - b = 0 -> b = a)
+      | [ H : ?x - ?y = 0, H' : ?x <> ?y |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a - b = 0 -> a = b)
+      | [ H : ?x * ?y * ?z <> 0, H' : ?y <> 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b c, a * b * c <> 0 -> a * c <> 0)
+      | [ H : ?x * ?y^3 <> 0, H' : ?y <> 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a * b <> 0 -> a <> 0)
+      | [ H : ?x * ?y^2 <> 0, H' : ?y <> 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a * b <> 0 -> a <> 0)
+      | [ H : ?x^2 - ?y^2 = ?z |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b c, a^2 - b^2 = c -> (a - b) * (a + b) = c)
+      | [ H : ?x + ?y * ?z = 0, H' : ?x <> Fopp (?z * ?y) |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b c, a + b * c = 0 -> a <> Fopp (c * b) -> False)
+      | [ H : ?x + ?x <> 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a, a + a <> 0 -> a <> 0)
+      | [ H : ?x - ?y <> 0, H' : ?y = ?x |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a - b <> 0 -> b <> a)
+      | [ H : ?x - ?y <> 0, H' : ?x = ?y |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a - b <> 0 -> a <> b)
+      | [ H : (?x + ?y)^2 - (?x^2 + ?y^2) = 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, (a + b)^2 - (a^2 + b^2) = 0 -> a * b = 0)
+      | [ H : (?x + ?y)^2 - (?x^2 + ?y^2) <> 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, (a + b)^2 - (a^2 + b^2) <> 0 -> a * b <> 0)
+      | [ H : ?x * ?y = 0, H' : ?x <> 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b, a * b = 0 -> a <> 0 -> b = 0)
+      | [ H : ?x * ?y = 0, H' : ?y <> 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b, a * b = 0 -> b <> 0 -> a = 0)
+      | [ H : ?x * ?y <> 0, H' : ?x <> 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a * b <> 0 -> b <> 0)
+      | [ H : ?x * ?y <> 0, H' : ?y <> 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a * b <> 0 -> a <> 0)
+      | [ H : ?x * ?y <> 0, H' : ?x = 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a * b <> 0 -> a <> 0)
+      | [ H : ?x * ?y <> 0, H' : ?y = 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, a * b <> 0 -> b <> 0)
+      | [ |- (?x + ?y)^2 - (?x^2 + ?y^2) = 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, a * b = 0 -> (a + b)^2 - (a^2 + b^2) = 0)
+      | [ |- (?x + ?y)^2 - (?x^2 + ?y^2) <> 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, a * b <> 0 -> (a + b)^2 - (a^2 + b^2) <> 0)
+      | [ |- (?x + ?y)^2 - ?x^2 - ?y^2 = 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, a * b = 0 -> (a + b)^2 - a^2 - b^2 = 0)
+      | [ |- (?x + ?y)^2 - ?x^2 - ?y^2 <> 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, a * b <> 0 -> (a + b)^2 - a^2 - b^2 <> 0)
+      | [ H : (?x + ?y)^2 - (?x^2 + ?y^2) = 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, (a + b)^2 - (a^2 + b^2) = 0 -> a * b = 0)
+      | [ H : (?x + ?y)^2 - (?x^2 + ?y^2) <> 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, (a + b)^2 - (a^2 + b^2) <> 0 -> a * b <> 0)
+      | [ H : (?x + ?y)^2 - ?x^2 - ?y^2 = 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, (a + b)^2 - a^2 - b^2 = 0 -> a * b = 0)
+      | [ H : (?x + ?y)^2 - ?x^2 - ?y^2 <> 0 |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b, (a + b)^2 - a^2 - b^2 <> 0 -> a * b <> 0)
+      | [ H : ?x = 0 |- ?x * ?y = 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, a = 0 -> a * b = 0)
+      | [ H : ?y = 0 |- ?x * ?y = 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, b = 0 -> a * b = 0)
+      | [ H : ?x <> 0 |- ?x * ?y <> 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, a <> 0 -> b <> 0 -> a * b <> 0)
+      | [ H : ?y <> 0 |- ?x * ?y <> 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, a <> 0 -> b <> 0 -> a * b <> 0)
+      | [ H : ?x <> 0 |- ?x * ?y = 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, b = 0 -> a * b = 0)
+      | [ H : ?y <> 0 |- ?x * ?y = 0 ]
+        => F.goal_apply_lemma_on pkg (forall a b, a = 0 -> a * b = 0)
+
+      | [ H : ?x - ?y * ?z <> ?w, H' : ?y = 1 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a - b * c <> d -> b = 1 -> a - c <> d) ltac:(fun H'' => constr:(fun Hv => H'' x y z w Hv H'))
+      | [ H : ?x - ?y * ?z = ?w, H' : ?y = 1 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a - b * c = d -> b = 1 -> a - c = d) ltac:(fun H'' => constr:(fun Hv => H'' x y z w Hv H'))
+      | [ H : ?x - ?y * ?z * ?w <> ?v, H' : ?y = 1 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d e, a - b * c * d <> e -> b = 1 -> a - c * d <> e) ltac:(fun H'' => constr:(fun Hv => H'' x y z w v Hv H'))
+      | [ H : ?x - ?y * ?z * ?w = ?v, H' : ?y = 1 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d e, a - b * c * d = e -> b = 1 -> a - c * d = e) ltac:(fun H'' => constr:(fun Hv => H'' x y z w v Hv H'))
+      | [ H : ?x - ?y * ?z^2 <> ?w, H' : ?z = 1 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a - b * c^2 <> d -> c = 1 -> a - b <> d) ltac:(fun H'' => constr:(fun Hv => H'' x y z w Hv H'))
+      | [ H : ?x - ?y * ?z^2 = ?w, H' : ?z = 1 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a - b * c^2 = d -> c = 1 -> a - b = d) ltac:(fun H'' => constr:(fun Hv => H'' x y z w Hv H'))
+
+      | [ H : ?x - ?y + ?z <> ?w, H' : ?y = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a - b + c <> d -> b = 0 -> a + c <> d) ltac:(fun H'' => constr:(fun Hv => H'' x y z w Hv H'))
+      | [ H : ?x - ?y <> ?z, H' : ?y = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a - b <> c -> b = 0 -> a <> c) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : ?x * ?y + ?z <> ?w, H' : ?x = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a * b + c <> d -> a = 0 -> c <> d) ltac:(fun H'' => constr:(fun Hv => H'' x y z w Hv H'))
+      | [ H : ?x * ?y + ?z <> ?w, H' : ?y = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a * b + c <> d -> b = 0 -> c <> d) ltac:(fun H'' => constr:(fun Hv => H'' x y z w Hv H'))
+      | [ H : ?x * ?y <> ?z, H' : ?x = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a * b <> c -> a = 0 -> c <> 0) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : ?x * ?y <> ?z, H' : ?y = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a * b <> c -> b = 0 -> c <> 0) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : Fopp (?x * ?y) <> ?z, H' : ?x = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, Fopp (a * b) <> c -> a = 0 -> c <> 0) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : Fopp (?x * ?y) <> ?z, H' : ?y = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, Fopp (a * b) <> c -> b = 0 -> c <> 0) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : ?x * ?y^3 = ?z, H' : ?y = 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b c, a * b^3 = c -> b = 0 -> c = 0)
+      | [ H : ?x * ?y = ?z, H' : ?x = 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b c, a * b = c -> a = 0 -> c = 0)
+      | [ H : ?x * ?y - ?z <> 0, H' : ?x = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a * b - c <> 0 -> a = 0 -> c <> 0) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : ?x * ?y - ?z <> 0, H' : ?y = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a * b - c <> 0 -> b = 0 -> c <> 0) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : ?x * ?y - ?z = 0, H' : ?x = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a * b - c = 0 -> a = 0 -> c = 0) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : ?x * ?y - ?z = 0, H' : ?y = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a * b - c = 0 -> b = 0 -> c = 0) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : ?x - ?y * ?z = 0, H' : ?z = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a - b * c = 0 -> c = 0 -> a = 0) ltac:(fun H'' => constr:(fun Hv => H'' x y z Hv H'))
+      | [ H : ?x * (?y * ?y^2) - ?z <> ?w |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b c d, a * (b * b^2) - c <> d -> a * (b^3) - c <> d)
+      | [ H : ?x * (?y * ?y^2) - ?z = ?w |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b c d, a * (b * b^2) - c = d -> a * (b^3) - c = d)
+      | [ H : ?x - (?y * ?y^2) * ?z <> ?w |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b c d, a - (b * b^2) * c <> d -> a - (b^3) * c <> d)
+      | [ H : ?x - (?y * ?y^2) * ?z = ?w |- _ ]
+        => F.apply_lemma_in_on pkg H (forall a b c d, a - (b * b^2) * c = d -> a - (b^3) * c = d)
+      | [ H : ?x * (?y * ?y^2) = 0, H' : ?y <> 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b, a * (b * b^2) = 0 -> b <> 0 -> a = 0)
+      | [ H : ?x * (?y * ?z) = 0, H' : ?z <> 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, c <> 0 -> a * (b * c) = 0 -> a * b = 0) ltac:(fun lem => constr:(lem x y z H'))
+      | [ H : ?x = ?y * ?z, H' : ?y = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a = b * c -> b = 0 -> a = 0) ltac:(fun lem => constr:(fun H => lem x y z H H'))
+      | [ H : ?x * ?y + ?z = 0, H' : ?x = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c, a = 0 -> a * b + c = 0 -> c = 0) ltac:(fun lem => constr:(lem x y z H'))
+      | [ H : ?x * (?y * ?y^2) - ?z * ?z^2 * ?w = 0, H' : ?x * ?y^3 + ?w * ?z^3 = 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a * b^3 + d * c^3 = 0 -> a * (b * b^2) - c * c^2 * d = 0 -> a * b^3 = 0) ltac:(fun lem => constr:(lem x y z w H'))
+      | [ H : ?x * Finv (?y^2) <> ?z * Finv (?w^2), Hy : ?y <> 0, Hw : ?w <> 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a * Finv (b^2) <> c * Finv (d^2) -> b <> 0 -> d <> 0 -> a * (d^2) <> c * (b^2)) ltac:(fun Hv => constr:(fun H => Hv x y z w H Hy Hw))
+      | [ H : ?x * Finv (?y^2) = ?z * Finv (?w^2), Hy : ?y <> 0, Hw : ?w <> 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a * Finv (b^2) = c * Finv (d^2) -> b <> 0 -> d <> 0 -> a * (d^2) = c * (b^2)) ltac:(fun Hv => constr:(fun H => Hv x y z w H Hy Hw))
+      | [ H : ?x * Finv (?y^3) = Fopp (?z * Finv (?w^3)), Hy : ?y <> 0, Hw : ?w <> 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a * Finv (b^3) = Fopp (c * Finv (d^3)) -> b <> 0 -> d <> 0 -> a * (d^3) = Fopp (c * (b^3))) ltac:(fun Hv => constr:(fun H => Hv x y z w H Hy Hw))
+      | [ H : ?x * Finv (?y^3) <> Fopp (?z * Finv (?w^3)), Hy : ?y <> 0, Hw : ?w <> 0 |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall a b c d, a * Finv (b^3) <> Fopp (c * Finv (d^3)) -> b <> 0 -> d <> 0 -> a * (d^3) <> Fopp (c * (b^3))) ltac:(fun Hv => constr:(fun H => Hv x y z w H Hy Hw))
+      | [ H : ?x = Fopp (?y * ?z), H' : ?x - ?z * ?y = 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b c, a = Fopp (b * c) -> a - c * b = 0 -> a = 0)
+      | [ H : ?x = Fopp (?y * ?z), H' : ?x - ?z * ?y = 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b c, a = Fopp (b * c) -> a - c * b = 0 -> a = 0)
+      | [ H : ?x <> Fopp (?y * ?z), H' : ?x - ?z * ?y = 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b c, a <> Fopp (b * c) -> a - c * b = 0 -> a * b * c <> 0)
+      | [ H : ?x <> Fopp (?y * ?z), H' : ?x - ?z * ?y = 0 |- _ ]
+        => F.apply2_lemma_in_on pkg H H' (forall a b c, a <> Fopp (b * c) -> a - c * b = 0 -> a * b * c <> 0)
+      | [ H : ?y^2 = ?c^3 + ?a * ?c * (?x^2)^2 + ?b * (?x^3)^2,
+              H' : ?y'^2 = ?c'^3 + ?a * ?c' * (?x'^2)^2 + ?b * (?x'^3)^2,
+                   H0 : ?c' * ?x^2 - ?c * ?x'^2 = 0
+          |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall X Y X' Y' A B C C', Y^2 = C^3 + A * C * (X^2)^2 + B * (X^3)^2 -> Y'^2 = C'^3 + A * C' * (X'^2)^2 + B * (X'^3)^2 -> C' * X^2 - C * X'^2 = 0 -> (Y' * (X^3))^2 - ((X'^3) * Y)^2 = 0) ltac:(fun Hv => constr:(fun H => Hv x y x' y' a b c c' H H' H0))
+      | [ H : ?y^2 = ?c^3 + ?a * ?c * (?x^2)^2 + ?b * (?x^3)^2,
+              H' : ?y'^2 = ?c'^3 + ?a * ?c' * (?x'^2)^2 + ?b * (?x'^3)^2,
+                   H0 : ?c' * ?x^2 = ?c * ?x'^2
+          |- _ ]
+        => F.apply_lemma_in_on' pkg H (forall X Y X' Y' A B C C', Y^2 = C^3 + A * C * (X^2)^2 + B * (X^3)^2 -> Y'^2 = C'^3 + A * C' * (X'^2)^2 + B * (X'^3)^2 -> C' * X^2 = C * X'^2 -> (Y' * (X^3))^2 - ((X'^3) * Y)^2 = 0) ltac:(fun Hv => constr:(fun H => Hv x y x' y' a b c c' H H' H0))
+
+      | [ H0 : ?x * ?y^3 = ?A * ?B^3,
+               H1 : ?x * ?z^3 - ?B^3 * ?D = 0,
+                    H2 : ?D * ?C^3 = ?w * ?z^3,
+                         H3 : ?A * ?C^3 - ?y^3 * ?w <> 0,
+                              Hz : ?z <> 0,
+                                   HB : ?B <> 0
+          |- _ ]
+        => exfalso; revert H0 H1 H2 H3 Hz HB;
+           generalize x, y, z, w, A, B, C, D;
+           F.goal_exact_lemma_on pkg
+      | [ H0 : ?x * ?y^3 = ?A * ?B^3,
+               H1 : ?x * ?z^3 - ?B^3 * ?D <> 0,
+                    H2 : ?D * ?C^3 = ?w * ?z^3,
+                         H3 : ?A * ?C^3 - ?y^3 * ?w = 0,
+                              Hy : ?y <> 0,
+                                   HC : ?C <> 0
+          |- _ ]
+        => exfalso; revert H0 H1 H2 H3 Hy HC;
+           generalize x, y, z, w, A, B, C, D;
+           F.goal_exact_lemma_on pkg
+      | [ H0 : ?x * ?y^2 = ?A * ?B^2,
+               H1 : ?x * ?z^2 - ?D * ?B^2 = 0,
+                    H2 : ?D * ?C^2 = ?w * ?z^2,
+                         H3 : ?A * ?C^2 - ?w * ?y^2 <> 0,
+                              Hz : ?z <> 0,
+                                   HB : ?B <> 0
+          |- _ ]
+        => exfalso; revert H0 H1 H2 H3 Hz HB;
+           generalize x, y, z, w, A, B, C, D;
+           F.goal_exact_lemma_on pkg
+      | [ H0 : ?A * ?B^2 = ?x * ?y^2,
+               H1 : ?x * ?z^2 - ?D * ?B^2 = 0,
+                    H2 : ?w * ?z^2 = ?D * ?C^2,
+                         H3 : ?A * ?C^2 - ?w * ?y^2 <> 0,
+                              Hz : ?z <> 0,
+                                   HB : ?B <> 0
+          |- _ ]
+        => exfalso; revert H0 H1 H2 H3 Hz HB;
+           generalize x, y, z, w, A, B, C, D;
+           F.goal_exact_lemma_on pkg
+      | [ H : ?x <> 0 |- context[Finv (?x^2)] ]
+        => F.with_lemma_on pkg (forall a, a <> 0 -> Finv (a^2) = (Finv a)^2) ltac:(fun H' => rewrite !(H' x H))
+      | [ H : ?x <> 0 |- context[Finv (?x^3)] ]
+        => F.with_lemma_on pkg (forall a, a <> 0 -> Finv (a^3) = (Finv a)^3) ltac:(fun H' => rewrite !(H' x H))
+      | [ H : ?x <> 0 |- context[Finv (Finv ?x)] ]
+        => F.with_lemma_on pkg (forall a, a <> 0 -> Finv (Finv a) = a) ltac:(fun H' => rewrite !(H' x H))
+      | [ |- context[Finv ((?x + ?y)^2 - ?x^2 - ?y^2)] ]
+        => F.with_lemma_on pkg (forall a b, (a + b)^2 - a^2 - b^2 = (1+1) * a * b) ltac:(fun H' => rewrite !(H' x y))
+      | [ |- context[Finv ((?x + ?y)^2 - (?x^2 + ?y^2))] ]
+        => F.with_lemma_on pkg (forall a b, (a + b)^2 - (a^2 + b^2) = (1+1) * a * b) ltac:(fun H' => rewrite !(H' x y))
+      | [ |- context[Finv (?x * ?y)] ]
+        => F.with_lemma_on
+             pkg
+             (forall a b, a * b <> 0 -> Finv (a * b) = Finv a * Finv b)
+             ltac:(fun H' => rewrite !(H' x y) by (clear_eq; fsatz))
+      | _ => idtac
+      end.
+    Local Ltac speed_up_fsatz :=
+      let fld := guess_field in
+      let pkg := F.get_package_on fld in
+      divisions_to_inverses fld;
+      repeat speed_up_fsatz_step_on pkg.
+    Local Ltac speed_up_fsatz_check :=
+      let fld := guess_field in
+      let pkg := F.get_package_on fld in
+      speed_up_fsatz_step_on pkg.
+
+    Local Ltac pre_pre_faster_t :=
+      prept; try contradiction; speed_up_fsatz.
+    Local Ltac pre_faster_t :=
+      pre_pre_faster_t; speed_up_fsatz_check; clean_up_speed_up_fsatz.
+    Local Ltac faster_t_noclear :=
+      pre_faster_t; fsatz.
+    Local Ltac faster_t :=
+      pre_faster_t;
+      try solve [ lazymatch goal with
+                  | [ |- _ = _ ] => clear_neq
+                  | _ => idtac
+                  end;
+                  fsatz ].
+
+    Hint Unfold Jacobian.double negb andb : points_as_coordinates.
+    Hint Unfold W.eq W.add Jacobian.to_affine Jacobian.of_affine Jacobian.add Jacobian.add_impl Jacobian.opp : points_as_coordinates.
+
+    (* These could go into Jacobian.v *)
+    Global Instance Proper_opp : Proper (eq ==> eq) opp. Proof. faster_t_noclear. Qed.
+
+    Lemma of_affine_add P Q
+      : eq (of_affine (W.add P Q)) (add (of_affine P) (of_affine Q)).
+    Proof. t. Qed.
+
+    Lemma add_opp (P : point) :
+      z_of (add P (opp P)) = 0.
+    Proof. faster_t_noclear. Qed.
+
+    Lemma add_comm (P Q : point) :
+      eq (add P Q) (add Q P).
+    Proof. faster_t_noclear. Qed.
+
+    Lemma add_zero_l (P Q : point) (H : z_of P = 0) :
+      eq (add P Q) Q.
+    Proof. faster_t. Qed.
+
+    Lemma add_zero_r (P Q : point) (H : z_of Q = 0) :
+      eq (add P Q) P.
+    Proof. faster_t. Qed.
+
+    Lemma add_double (P Q : point) :
+      eq P Q ->
+      eq (add P Q) (double P).
+    Proof. faster_t_noclear. Qed.
+
+    (* This uses assumptions not present in Jacobian.v,
+       namely char_ge_12 and discriminant_nonzero. *)
+    Lemma add_assoc (P Q R : point) :
+      eq (add (add P Q) R) (add P (add Q R)).
+    Proof.
+      rewrite <- (Jacobian.of_affine_to_affine (add (add _ _) _)).
+      rewrite <- (Jacobian.of_affine_to_affine (add _ (add _ _))).
+      rewrite (Jacobian.to_affine_add _ R).
+      rewrite (Jacobian.to_affine_add P (add Q R)).
+      do 2 rewrite of_affine_add.
+      rewrite (Jacobian.to_affine_add P Q).
+      rewrite (Jacobian.to_affine_add Q R).
+      do 2 rewrite <- of_affine_add.
+      destruct (@W.commutative_group F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec char_ge_3 char_ge_12 ltac:(cbv [id] in *; fsatz)).
+      destruct commutative_group_group. destruct group_monoid.
+      destruct monoid_is_associative.
+      rewrite associative.
+      reflexivity.
+    Qed.
+
+    Definition is_point (P : F*F*F) : Prop :=
+      let '(X,Y,Z) := P in
+      if dec (Z=0) then True else Y^2 = X^3 + a*X*(Z^2)^2 + b*(Z^3)^2.
+
+    Definition co_z (P Q : point) : Prop :=
+      z_of P = z_of Q.
+
+    Lemma opp_co_z (P : point) :
+      co_z P (opp P).
+    Proof. unfold co_z; faster_t. Qed.
+
+    Program Definition make_co_z (P Q : point) (HQaff : z_of Q = 1) : point*point :=
+      match proj1_sig P, proj1_sig Q return (F*F*F)*(F*F*F) with
+      | (x1, y1, z1), (x2, y2, z2) =>
+          let t1 := x2 in
+          let t2 := y2 in
+          let t3 := z1 in
+          let t4 := t3^2 in
+          let t1 := t4 * t1 in
+          let t2 := t4 * t2 in
+          let t2 := t3 * t2 in
+          (P, (t1, t2, t3))
+      end.
+    Next Obligation. Proof. faster_t. Qed.
+    Next Obligation. Proof. faster_t. Qed.
+
+    Hint Unfold is_point co_z make_co_z : points_as_coordinates.
+
+    Lemma make_co_z_correct (P Q : point) (HQaff : z_of Q = 1)
+      (HPnz : z_of P <> 0) :
+      let '(A, B) := make_co_z P Q HQaff in
+      eq P A
+      /\ eq Q B
+      /\ co_z A B.
+    Proof. faster_t. Qed.
+
+    (* Scalar Multiplication on Weierstraß Elliptic Curves from Co-Z Arithmetic *)
+    (* Goundar, Joye, Miyaji, Rivain, Vanelli *)
+    (* Algorithm 19 Co-Z addition with update (register allocation) *)
+    (* 6 field registers, 5M, 2S *)
+    Program Definition zaddu (P Q : point) (H : co_z P Q) : point * point :=
+      match proj1_sig P, proj1_sig Q return (F*F*F)*(F*F*F) with
+      | (x1, y1, z1), (x2, y2, z2) =>
+          let t1 := x1 in
+          let t2 := y1 in
+          let t3 := z1 in
+          let t4 := x2 in
+          let t5 := y2 in
+          let t6 := t1 - t4 in
+          let t3 := t3 * t6 in
+          let t6 := t6^2 in
+          let t1 := t1 * t6 in
+          let t6 := t6 * t4 in
+          let t5 := t2 - t5 in
+          let t4 := t5^2 in
+          let t4 := t4 - t1 in
+          let t4 := t4 - t6 in
+          let t6 := t1 - t6 in
+          let t2 := t2 * t6 in
+          let t6 := t1 - t4 in
+          let t5 := t5 * t6 in
+          let t5 := t5 - t2 in
+          ((t4, t5, t3), (t1, t2, t3))
+      end.
+    Next Obligation. Proof. faster_t_noclear. Qed.
+    Next Obligation. Proof. faster_t. Qed.
+
+    Hint Unfold zaddu : points_as_coordinates.
+
+    (* ZADDU(P, Q) = (P + Q, P) if P <> Q, Q <> -P *)
+    Lemma zaddu_correct (P Q : point) (H : co_z P Q)
+      (Hneq : x_of P <> x_of Q):
+      let '(R1, R2) := zaddu P Q H in
+      eq (add P Q) R1 /\ eq P R2 /\ co_z R1 R2.
+    Proof. faster_t_noclear. Qed.
+
+    Lemma zaddu_correct_alt (P Q : point) (H : co_z P Q) :
+      let '(R1, R2) := zaddu P Q H in
+      z_of R1 <> 0 ->
+      eq (add P Q) R1 /\ eq P R2 /\ co_z R1 R2.
+    Proof.
+      generalize (zaddu_correct P Q H).
+      rewrite (surjective_pairing (zaddu _ _ _)).
+      intros A B. apply A.
+      clear -B. t.
+    Qed.
+
+    Lemma zaddu_correct0 (P : point) :
+      let '(R1, R2) := zaddu P (opp P) (opp_co_z P) in
+      z_of R1 = 0 /\ co_z R1 R2.
+    Proof. faster_t_noclear. Qed.
+
+    (* Scalar Multiplication on Weierstraß Elliptic Curves from Co-Z Arithmetic *)
+    (* Goundar, Joye, Miyaji, Rivain, Vanelli *)
+    (* Algorithm 20 Conjugate co-Z addition (register allocation) *)
+    (* 7 field registers, 6M, 3S*)
+    Program Definition zaddc (P Q : point) (H : co_z P Q) : point * point :=
+      match proj1_sig P, proj1_sig Q return (F*F*F)*(F*F*F) with
+      | (x1, y1, z1), (x2, y2, z2) =>
+          let t1 := x1 in
+          let t2 := y1 in
+          let t3 := z1 in
+          let t4 := x2 in
+          let t5 := y2 in
+          let t6 := t1 - t4 in
+          let t3 := t3 * t6 in
+          let t6 := t6^2 in
+          let t7 := t1 * t6 in
+          let t6 := t6 * t4 in
+          let t1 := t2 + t5 in
+          let t4 := t1^2 in
+          let t4 := t4 - t7 in
+          let t4 := t4 - t6 in
+          let t1 := t2 - t5 in
+          let t1 := t1^2 in
+          let t1 := t1 - t7 in
+          let t1 := t1 - t6 in
+          let t6 := t6 - t7 in
+          let t6 := t6 * t2 in
+          let t2 := t2 - t5 in
+          let t5 := t5 + t5 in
+          let t5 := t2 + t5 in
+          let t7 := t7 - t4 in
+          let t5 := t5 * t7 in
+          let t5 := t5 + t6 in
+          let t7 := t4 + t7 in
+          let t7 := t7 - t1 in
+          let t2 := t2 * t7 in
+          let t2 := t2 + t6 in
+          ((t1, t2, t3), (t4, t5, t3))
+      end.
+    Next Obligation. Proof. faster_t_noclear. Qed.
+    Next Obligation. Proof. faster_t_noclear. Qed.
+
+    Hint Unfold zaddc : points_as_coordinates.
+    (* ZADDC(P, Q) = (P + Q, P - Q) if P <> Q, Q <> -P *)
+    Lemma zaddc_correct (P Q : point) (H : co_z P Q)
+      (Hneq : x_of P <> x_of Q):
+      let '(R1, R2) := zaddc P Q H in
+      eq (add P Q) R1 /\ eq (add P (opp Q)) R2 /\ co_z R1 R2.
+    Proof. faster_t_noclear. Qed.
+
+    Lemma zaddc_correct_alt (P Q : point) (H : co_z P Q) :
+      let '(R1, R2) := zaddc P Q H in
+      z_of R1 <> 0 ->
+      eq (add P Q) R1 /\ eq (add P (opp Q)) R2 /\ co_z R1 R2.
+    Proof.
+      generalize (zaddc_correct P Q H).
+      rewrite (surjective_pairing (zaddc _ _ _)).
+      intros A B. apply A.
+      clear -B. t.
+    Qed.
+
+    Lemma zaddc_correct0 (P : point) :
+      let '(R1, R2) := zaddc P (opp P) (opp_co_z P) in
+      z_of R1 = 0 /\ co_z R1 R2.
+    Proof. faster_t_noclear. Qed.
+
+    (* Scalar Multiplication on Weierstraß Elliptic Curves from Co-Z Arithmetic *)
+    (* Goundar, Joye, Miyaji, Rivain, Vanelli *)
+    (* Algorithm 21 Co-Z doubling with update (register allocation) *)
+    (* 6 field registers, 1M, 5S *)
+    Program Definition dblu (P : point) (H: z_of P = 1) : point * point :=
+      match proj1_sig P return (F*F*F)*(F*F*F) with
+      | (x1, y1, _) =>
+          let t0 :=  a in
+          let t1 := x1 in
+          let t2 := y1 in
+          let t3 := t2 + t2 in
+          let t2 := t2^2 in
+          let t4 := t1 + t2 in
+          let t4 := t4^2 in
+          let t5 := t1^2 in
+          let t4 := t4 - t5 in
+          let t2 := t2^2 in
+          let t4 := t4 - t2 in
+          let t1 := t4 + t4 in
+          let t0 := t0 + t5 in
+          let t5 := t5 + t5 in
+          let t0 := t0 + t5 in
+          let t4 := t0^2 in
+          let t5 := t1 + t1 in
+          let t4 := t4 - t5 in
+          let t2 := 8 * t2 in
+          let t5 := t1 - t4 in
+          let t5 := t5 * t0 in
+          let t5 := t5 - t2 in
+          ((t4, t5, t3), (t1, t2, t3))
+      end.
+    Next Obligation. Proof. faster_t. Qed.
+    Next Obligation. Proof. faster_t. Qed.
+
+    Hint Unfold dblu : points_as_coordinates.
+
+    (* DBLU(P) = (2P, P) when Z(P) = 1 *)
+    Lemma dblu_correct (P : point) (H : z_of P = 1)
+      (Hynz : y_of P <> 0) :
+      let '(R1, R2) := dblu P H in
+      eq (double P) R1 /\ eq P R2 /\ co_z R1 R2.
+    Proof. faster_t. Qed.
+
+    Lemma dblu_correct_alt (P : point) (H : z_of P = 1) :
+      let '(R1, R2) := dblu P H in
+      z_of R1 <> 0 ->
+      eq (double P) R1 /\ eq P R2 /\ co_z R1 R2.
+    Proof.
+      generalize (dblu_correct P H).
+      rewrite (surjective_pairing (dblu _ _)).
+      intros A B. apply A.
+      clear -B. t.
+    Qed.
+
+    Lemma dblu_correct0 (P : point) (H : z_of P = 1)
+      (Hyz : y_of P = 0) :
+      let '(R1, R2) := dblu P H in
+      z_of R1 = 0 /\ co_z R1 R2.
+    Proof. faster_t. Qed.
+
+    (* Scalar Multiplication on Weierstraß Elliptic Curves from Co-Z Arithmetic *)
+    (* Goundar, Joye, Miyaji, Rivain, Vanelli *)
+    (* Algorithm 22 Co-Z tripling with update (register allocation) *)
+    (* 6M, 7S *)
+    Program Definition tplu (P : point) (H : z_of P = 1) : point * point :=
+      zaddu (snd (dblu P H)) (fst (dblu P H)) _.
+    Next Obligation. faster_t. Qed.
+
+    Hint Unfold tplu : points_as_coordinates.
+
+    Lemma tplu_correct0 (P : point) (H : z_of P = 1) (Hyz : y_of P = 0) :
+      let '(R1, R2) := tplu P H in
+      z_of R1 = 0 /\ co_z R1 R2.
+    Proof. faster_t. Qed.
+
+    Lemma tplu_correct (P : point) (H : z_of P = 1) (Hynz : y_of P <> 0) :
+      let '(R1, R2) := tplu P H in
+      z_of R1 <> 0 ->
+      eq (add (double P) P) R1 /\ eq P R2 /\ co_z R1 R2.
+    Proof.
+      rewrite (surjective_pairing (tplu P H)).
+      unfold co_z, tplu.
+      generalize (zaddu_correct_alt (snd (dblu P H)) (fst (dblu P H)) (tplu_obligation_1 P H)).
+      destruct (zaddu (snd (dblu P H)) (fst (dblu P H)) (tplu_obligation_1 P H)) as [R1 R2] eqn:Hzaddu.
+      intros A B. specialize (A B). destruct A as [A1 [A2 A3] ].
+      generalize (dblu_correct P H Hynz).
+      rewrite (surjective_pairing (dblu P H)). intros [B1 [B2 B3] ].
+      simpl; repeat split; [| |assumption].
+      - rewrite B1. rewrite <- A1.
+        rewrite B2 at 2. apply add_comm.
+      - etransitivity; eauto.
+    Qed.
+
+    (* Scalar Multiplication on Weierstraß Elliptic Curves from Co-Z Arithmetic *)
+    (* Goundar, Joye, Miyaji, Rivain, Vanelli *)
+    (* Co-Z doubling-addition with update (naive version) *)
+    (* 11M, 5S *)
+    Program Definition zdau_naive (P Q : point) (H : co_z P Q) :=
+      zaddc (fst (zaddu P Q H)) (snd (zaddu P Q H)) _.
+    Next Obligation. Proof. t. Qed.
+
+    (* ZDAU(P, Q) = (2P + Q, Q) if P <> Q, Q <> -P *)
+    Lemma zdau_naive_correct (P Q : point) (H : co_z P Q)
+      (Hneq : x_of P <> x_of Q) :
+      let '(R1, R2) := zdau_naive P Q H in
+      z_of R1 <> 0 ->
+      eq (add (double P) Q) R1 /\ eq Q R2 /\ co_z R1 R2.
+    Proof.
+      destruct (zdau_naive P Q H) as [R1 R2] eqn:HR.
+      intros HR1. unfold zdau_naive in HR.
+      generalize (zaddc_correct_alt (fst (zaddu P Q H)) (snd (zaddu P Q H)) (zdau_naive_obligation_1 P Q H)). rewrite HR.
+      intros A. specialize (A HR1).
+      destruct A as (A1 & A2 & A3).
+      generalize (zaddu_correct P Q H Hneq).
+      rewrite (surjective_pairing (zaddu P Q H)).
+      intros (B1 & B2 & B3).
+      repeat split; auto.
+      - rewrite <- add_double, <- A1, <- B1, <- B2; [|reflexivity].
+        rewrite add_assoc, add_comm. reflexivity.
+      - rewrite <- A2, <- B1, <- B2.
+        rewrite (add_comm P Q).
+        rewrite add_assoc. rewrite add_zero_r; [reflexivity|apply add_opp].
+    Qed.
+
+    Lemma zdau_naive_correct_alt (P Q : point) (H : co_z P Q)
+      (Hneq : x_of P <> x_of Q)
+      (Hneq2: x_of (fst (zaddu P Q H)) <> x_of (snd (zaddu P Q H))) :
+      let '(R1, R2) := zdau_naive P Q H in
+      eq (add (double P) Q) R1 /\ eq Q R2 /\ co_z R1 R2.
+    Proof.
+      destruct (zdau_naive P Q H) as [R1 R2] eqn:HR.
+      unfold zdau_naive in HR.
+      generalize (zaddc_correct (fst (zaddu P Q H)) (snd (zaddu P Q H)) (zdau_naive_obligation_1 P Q H) Hneq2). rewrite HR.
+      intros (A1 & A2 & A3).
+      generalize (zaddu_correct P Q H Hneq).
+      rewrite (surjective_pairing (zaddu P Q H)).
+      intros (B1 & B2 & B3).
+      repeat split; auto.
+      - rewrite <- add_double, <- A1, <- B1, <- B2; [|reflexivity].
+        rewrite add_assoc, add_comm. reflexivity.
+      - rewrite <- A2, <- B1, <- B2.
+        rewrite (add_comm P Q).
+        rewrite add_assoc. rewrite add_zero_r; [reflexivity|apply add_opp].
+    Qed.
+
+    (* Scalar Multiplication on Weierstraß Elliptic Curves from Co-Z Arithmetic *)
+    (* Goundar, Joye, Miyaji, Rivain, Vanelli *)
+    (* Algorithm 23 Co-Z doubling-addition with update (register allocation) *)
+    (* 8 field registers, 9M, 7S (instead of 11M, 5S) *)
+    Program Definition zdau (P Q : point) (H : co_z P Q) : point * point :=
+      match proj1_sig P, proj1_sig Q return (F*F*F)*(F*F*F) with
+      | (x1, y1, z1), (x2, y2, z2) =>
+          let t1 := x1 in
+          let t2 := y1 in
+          let t3 := z1 in
+          let t4 := x2 in
+          let t5 := y2 in
+          let t6 := t1 - t4 in
+          let t7 := t6^2 in
+          let t1 := t1 * t7 in
+          let t4 := t4 * t7 in
+          let t5 := t2 - t5 in
+          let t8 := t1 - t4 in
+          let t2 := t2 * t8 in
+          let t2 := t2 + t2 in
+          let t8 := t5^2 in
+          let t4 := t8 - t4 in
+          let t4 := t4 - t1 in
+          let t4 := t4 - t1 in
+          let t6 := t4 + t6 in
+          let t6 := t6^2 in
+          let t6 := t6 - t7 in
+          let t5 := t5 - t4 in
+          let t5 := t5^2 in
+          let t5 := t5 - t8 in
+          let t5 := t5 - t2 in
+          let t7 := t4^2 in
+          let t5 := t5 - t7 in
+          let t8 := t7 + t7 in
+          let t8 := t8 + t8 in
+          let t6 := t6 - t7 in
+          let t3 := t3 * t6 in
+          let t6 := t1 * t8 in
+          let t1 := t1 + t4 in
+          let t8 := t8 * t1 in
+          let t7 := t2 + t5 in
+          let t2 := t5 - t2 in
+          let t1 := t8 - t6 in
+          let t5 := t5 * t1 in
+          let t6 := t6 + t8 in
+          let t1 := t2^2 in
+          let t1 := t1 - t6 in
+          let t4 := t8 - t1 in
+          let t2 := t2 * t4 in
+          let t2 := t2 - t5 in
+          let t4 := t7^2 in
+          let t4 := t4 - t6 in
+          let t8 := t8 - t4 in
+          let t7 := t7 * t8 in
+          let t5 := t7 - t5 in
+          ((t1, t2, t3), (t4, t5, t3))
+      end.
+    Next Obligation. Proof. faster_t_noclear. Qed.
+    Next Obligation. Proof. faster_t_noclear. Qed.
+
+    Hint Unfold zdau : points_as_coordinates.
+
+    Lemma zdau_naive_eq_zdau (P Q : point) (H : co_z P Q) :
+      let '(R1, R2) := zdau_naive P Q H in
+      let '(S1, S2) := zdau P Q H in
+      eq R1 S1 /\ eq R2 S2.
+    Proof. faster_t. all: try fsatz. Qed.
+
+    (* Direct proof is intensive, which is why we need the naive implementation *)
+    Lemma zdau_correct (P Q : point) (H : co_z P Q)
+      (Hneq : x_of P <> x_of Q) :
+      let '(R1, R2) := zdau P Q H in
+      z_of R1 <> 0 ->
+      eq (add (double P) Q) R1 /\ eq Q R2 /\ co_z R1 R2.
+    Proof.
+      generalize (zdau_naive_correct P Q H Hneq).
+      generalize (zdau_naive_eq_zdau P Q H).
+      rewrite (surjective_pairing (zdau_naive _ _ _)).
+      rewrite (surjective_pairing (zdau _ _ _)).
+      intros [A1 A2] X Y.
+      assert (Z: forall A B, eq A B -> z_of A <> 0 -> z_of B <> 0) by (clear; intros; faster_t).
+      repeat split.
+      - rewrite <- A1. apply X. eapply Z; eauto. symmetry; assumption.
+      - rewrite <- A2. apply X. eapply Z; eauto. symmetry; assumption.
+      - clear -H. t.
+    Qed.
+
+    Lemma zdau_correct_alt (P Q : point) (H : co_z P Q)
+      (Hneq : x_of P <> x_of Q)
+      (Hneq2: x_of (fst (zaddu P Q H)) <> x_of (snd (zaddu P Q H))) :
+      let '(R1, R2) := zdau P Q H in
+      eq (add (double P) Q) R1 /\ eq Q R2 /\ co_z R1 R2.
+    Proof.
+      generalize (zdau_naive_correct_alt P Q H Hneq Hneq2).
+      generalize (zdau_naive_eq_zdau P Q H).
+      rewrite (surjective_pairing (zdau_naive _ _ _)).
+      rewrite (surjective_pairing (zdau _ _ _)).
+      intros [A1 A2] [A3 [A4 A5] ].
+      assert (Z: forall A B, eq A B -> z_of A <> 0 -> z_of B <> 0) by (clear; intros; faster_t).
+      repeat split.
+      - rewrite <- A1. auto.
+      - rewrite <- A2. auto.
+      - clear -H. t.
+    Qed.
+
+    Lemma zdau_correct0 (P : point) :
+      let '(R1, R2) := zdau P (opp P) (opp_co_z P) in
+      z_of R1 = 0 /\ co_z R1 R2.
+    Proof. faster_t. Qed.
+  End Co_Z.
+End Jacobian.
diff --git a/src/Curves/Weierstrass/Jacobian.v b/src/Curves/Weierstrass/Jacobian/Jacobian.v
similarity index 99%
rename from src/Curves/Weierstrass/Jacobian.v
rename to src/Curves/Weierstrass/Jacobian/Jacobian.v
index fd3b128b6e..562ee287db 100644
--- a/src/Curves/Weierstrass/Jacobian.v
+++ b/src/Curves/Weierstrass/Jacobian/Jacobian.v
@@ -539,7 +539,7 @@ Module Jacobian.
       eq (add P Q) (add_mixed P Q H).
     Proof. faster_t. Qed.
 
-    Lemma Proper_add : Proper (eq ==> eq ==> eq) add. Proof. faster_t_noclear. Qed.
+    Global Instance Proper_add : Proper (eq ==> eq ==> eq) add. Proof. faster_t_noclear. Qed.
     Import BinPos.
     Lemma to_affine_add P Q
       : W.eq (to_affine (add P Q)) (W.add (to_affine P) (to_affine Q)).
diff --git a/src/Curves/Weierstrass/Jacobian/ScalarMult.v b/src/Curves/Weierstrass/Jacobian/ScalarMult.v
new file mode 100644
index 0000000000..2fe2dba939
--- /dev/null
+++ b/src/Curves/Weierstrass/Jacobian/ScalarMult.v
@@ -0,0 +1,810 @@
+Require Import Coq.Classes.Morphisms.
+Require Import Crypto.Spec.WeierstrassCurve Crypto.Algebra.ScalarMult.
+Require Import Crypto.Curves.Weierstrass.Jacobian.Jacobian.
+Require Import Crypto.Curves.Weierstrass.Affine Crypto.Curves.Weierstrass.AffineProofs.
+Require Import Crypto.Curves.Weierstrass.Jacobian.CoZ.
+Require Import Crypto.Util.Decidable Crypto.Algebra.Field.
+Require Import Crypto.Util.Tactics.BreakMatch.
+Require Import Crypto.Util.Tactics.DestructHead.
+Require Import Crypto.Util.Tactics.SpecializeBy.
+Require Import Crypto.Util.Tactics.SetoidSubst.
+Require Import Crypto.Util.Notations Crypto.Util.LetIn.
+Require Import Crypto.Util.Sum Crypto.Util.Prod Crypto.Util.Sigma.
+Require Import Crypto.Util.FsatzAutoLemmas.
+Require Import Crypto.Util.Loops.
+Require Import Crypto.Util.ZUtil.Notations.
+Require Import Crypto.Util.ZUtil.Tactics.LtbToLt.
+Require Import Crypto.Util.ZUtil.Shift.
+Require Import Crypto.Util.ZUtil.Peano.
+Require Import Crypto.Util.Tuple.
+Require Import Coq.ZArith.ZArith.
+Require Import Coq.micromega.Lia.
+
+Module Z.
+(* Note: ideally we would contribute this to Coq *)
+Lemma mod_pow2_succ_r a n (Hi : Z.le 0 n) : (a mod 2 ^ Z.succ n = a mod 2 ^ n + Z.b2z (Z.testbit a n) * 2 ^ n)%Z.
+Proof.
+  clear -Hi.
+  rewrite Z_div_mod_eq_full with (a:=(a mod _)%Z) (b:=(2^n)%Z).
+  rewrite <-(Z.mod_small ((a mod 2 ^ Z.succ n / 2 ^ n)) 2); cycle 1.
+  { rewrite Z.pow_succ_r by lia; Z.div_mod_to_equations; nia. }
+  rewrite <-Z.mod_pow2_bits_low with (n:=Z.succ n) by lia.
+  rewrite <-Znumtheory.Zmod_div_mod; cycle 1.
+  { lia. } { lia. } { exists 2%Z. rewrite <-Z.pow_succ_r; f_equal; lia. }
+  rewrite Z.testbit_spec'; lia.
+Qed.
+End Z.
+
+Module ScalarMult.
+  Section JoyeDoubleAddDecomposition.
+    (* Joye's double-add ladder for computing [n]P basically tracks the
+       values of two sequences [SS i]P and [TT i]P.
+       This section proves some properties on the sequences SS and TT.
+     *)
+    Variables n : Z.
+    Local Coercion Z.of_nat : nat >-> Z.
+
+    Definition SS (i : nat) : Z := n mod 2 ^ i.
+    Definition TT (i : nat) : Z := 2^i - SS i.
+
+    Lemma SS_plus_TT k : (SS k + TT k = 2^(Z.of_nat k))%Z.
+    Proof. cbv [SS TT]; lia. Qed.
+
+    Lemma SS_succ i : SS (S i) = if (Z.testbit n (Z.of_nat i)) then (2 * SS i + TT i)%Z else SS i.
+    Proof.
+      cbv [TT SS].
+      pose proof Z.mod_pow2_succ_r n i; case Z.testbit in *; cbv [Z.b2z] in *; lia.
+    Qed.
+
+    Lemma TT_succ i : TT (S i) = if (Z.testbit n (Z.of_nat i)) then TT i else (2 * TT i + SS i)%Z.
+    Proof.
+      cbv [TT SS].
+      rewrite Nat2Z.inj_succ by lia.
+      pose proof Z.mod_pow2_succ_r n i; case Z.testbit in *; cbv [Z.b2z] in *; rewrite ?Z.pow_succ_r in *; lia.
+    Qed.
+
+    Lemma SS_monotone1 (k : nat) : (SS 1 <= SS (S k))%Z.
+    Proof.
+      cbv [SS].
+      rewrite Z.pow_1_r, Nat2Z.inj_succ, Z.pow_succ_r, Znumtheory.Zmod_div_mod with (n:=2) (m:=(2*2^k)%Z)
+        by (try exists (2^k)%Z; lia).
+      apply Z.mod_bound_pos_le; Z.div_mod_to_equations; lia.
+    Qed.
+
+    Lemma TT_monotone1 (k : nat) : (TT 1 <= TT (S k))%Z.
+    Proof.
+      cbv [TT SS]; repeat rewrite Nat2Z.inj_succ. set (z := Z.succ k).
+      enough (n mod 2 ^ z < 2 ^ z + (n mod 2 ^ Z.succ 0 - 1))%Z by lia; rewrite Z.pow_1_r.
+      enough (n mod 2 = 0 -> n mod 2 ^ z <> Z.pred (2^z))%Z by try (Z.div_mod_to_equations; nia); intros ? X.
+      apply (f_equal (fun x => Z.testbit x 0)) in X.
+      rewrite <-Z.ones_equiv, Z.ones_spec_low, Z.mod_pow2_bits_low, Z.bit0_eqb in X; lia.
+    Qed.
+
+    Lemma SSn scalarbitsz (Hn : (0 <= n < 2^scalarbitsz)%Z) (Hscalarbitsz : (0 < scalarbitsz)%Z) :
+      SS (Z.to_nat scalarbitsz) = n :> Z.
+    Proof.
+      cbv [SS]. rewrite Z2Nat.id by lia.
+      apply Z.mod_small; trivial.
+    Qed.
+  End JoyeDoubleAddDecomposition.
+
+  Section JoyeLadder.
+    Context {F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv} {a b:F}
+            {field:@Algebra.Hierarchy.field F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv}
+            {char_ge_3:@Ring.char_ge F Feq Fzero Fone Fopp Fadd Fsub Fmul (BinNat.N.succ_pos (BinNat.N.two))}
+            {Feq_dec:DecidableRel Feq}.
+    Local Infix "=" := Feq : type_scope. Local Notation "a <> b" := (not (a = b)) : type_scope.
+    Local Notation "0" := Fzero. Local Notation "1" := Fone.
+    Local Infix "+" := Fadd. Local Infix "-" := Fsub.
+    Local Infix "*" := Fmul. Local Infix "/" := Fdiv.
+    Local Notation "x ^ 2" := (x*x). Local Notation "x ^ 3" := (x^2*x).
+    Local Notation "2" := (1+1). Local Notation "4" := (2+1+1).
+    Local Notation "8" := (4+4). Local Notation "27" := (4*4 + 4+4 +1+1+1).
+    Local Notation "'∞'" := (@W.zero F Feq Fadd Fmul a b).
+    Local Notation Wpoint := (@W.point F Feq Fadd Fmul a b).
+    Context {char_ge_12:@Ring.char_ge F Feq Fzero Fone Fopp Fadd Fsub Fmul 12%positive}
+            {discriminant_nonzero:id(4*a*a*a + 27*b*b <> 0)}.
+    Local Notation Wopp := (@W.opp F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation Wadd := (@W.add F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv field Feq_dec char_ge_3 a b).
+    Local Notation Wzero := (@W.zero F Feq Fadd Fmul a b).
+    Local Notation Weq := (@W.eq F Feq Fadd Fmul a b).
+    Local Notation Wgroup := (@W.commutative_group F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec char_ge_3 char_ge_12 discriminant_nonzero).
+    Local Notation point := (@Jacobian.point F Feq Fzero Fadd Fmul a b Feq_dec).
+    Local Notation eq := (@Jacobian.eq F Feq Fzero Fadd Fmul a b Feq_dec).
+    Local Notation x_of := (@Jacobian.x_of F Feq Fzero Fadd Fmul a b Feq_dec).
+    Local Notation y_of := (@Jacobian.y_of F Feq Fzero Fadd Fmul a b Feq_dec).
+    Local Notation z_of := (@Jacobian.z_of F Feq Fzero Fadd Fmul a b Feq_dec).
+    Local Notation add := (@Jacobian.add F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field char_ge_3 Feq_dec).
+    Local Notation opp := (@Jacobian.opp F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation co_z := (@Jacobian.co_z F Feq Fzero Fadd Fmul a b Feq_dec).
+    Local Notation make_co_z := (@Jacobian.make_co_z F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation of_affine := (@Jacobian.of_affine F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation to_affine := (@Jacobian.to_affine F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation zero := (of_affine Wzero).
+    Local Notation dblu := (@Jacobian.dblu F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation tplu := (@Jacobian.tplu F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation zaddu := (@Jacobian.zaddu F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation zdau := (@Jacobian.zdau F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field Feq_dec).
+    Local Notation is_point := (@Jacobian.is_point F Feq Fzero Fadd Fmul a b Feq_dec).
+
+    Ltac prept_step_opt :=
+      match goal with
+      | [ H : ?T |- ?T ] => exact H
+      | [ |- ?x = ?x ] => reflexivity
+      | [ H : ?T, H' : ~?T |- _ ] => solve [ exfalso; apply H', H ]
+      end.
+
+    Ltac prept_step :=
+      match goal with
+      | _ => progress prept_step_opt
+      | _ => progress intros
+      (*| _ => progress specialize_by_assumption*)
+      (*| _ => progress specialize_by trivial*)
+      | _ => progress cbv [proj1_sig fst snd] in *
+      | _ => progress autounfold with points_as_coordinates in *
+      | _ => progress destruct_head'_True
+      | _ => progress destruct_head'_unit
+      | _ => progress destruct_head'_prod
+      | _ => progress destruct_head'_sig
+      | _ => progress destruct_head'_and
+      | _ => progress destruct_head'_sum
+      | _ => progress destruct_head'_bool
+      | _ => progress destruct_head'_or
+      | H: context[@dec ?P ?pf] |- _ => destruct (@dec P pf)
+      | |- context[@dec ?P ?pf]      => destruct (@dec P pf)
+      | |- ?P => lazymatch type of P with Prop => split end
+      end.
+    Ltac prept := repeat prept_step.
+    Ltac t := prept; trivial; try contradiction; fsatz.
+
+    Create HintDb points_as_coordinates discriminated.
+    Hint Unfold Proper respectful Reflexive Symmetric Transitive : points_as_coordinates.
+    Hint Unfold Jacobian.point Jacobian.eq W.eq W.point W.coordinates proj1_sig fst snd : points_as_coordinates.
+
+    Hint Unfold Jacobian.x_of Jacobian.y_of Jacobian.z_of Jacobian.opp Jacobian.co_z : points_as_coordinates.
+
+    Lemma co_z_comm (A B : point) (H : co_z A B) :
+      co_z B A.
+    Proof. t. Qed.
+
+    Definition co_z_points : Type := { '(A, B) | co_z A B }.
+
+    Program Definition make_co_z_points (P Q : point) (HQaff : z_of Q = 1) : co_z_points :=
+      make_co_z P Q HQaff.
+    Next Obligation. Proof. t. Qed.
+
+    Program Definition cswap_co_z_points (b : bool) (AB : co_z_points) : co_z_points :=
+      exist _
+        (let '(A, B) := proj1_sig AB
+        in if b then (B, A) else (A, B))
+        _.
+    Next Obligation. Proof. generalize (proj2_sig AB); rewrite (surjective_pairing (proj1_sig _)); destruct b0; [apply co_z_comm|auto]. Qed.
+
+    Program Definition zdau_co_z_points (AB : co_z_points) : co_z_points :=
+      exist _ (let '(A, B) := proj1_sig AB in zdau A B _) _.
+    Next Obligation. Proof. generalize (proj2_sig AB). rewrite <- Heq_anonymous. auto. Qed.
+    Next Obligation. Proof. destruct AB as ((A & B) & HAB). simpl. t. Qed.
+
+    Program Definition zaddu_co_z_points (AB : co_z_points) : co_z_points :=
+      exist _ (let '(A, B) := proj1_sig AB in zaddu A B _) _.
+    Next Obligation. Proof. generalize (proj2_sig AB); rewrite <- Heq_anonymous. auto. Qed.
+    Next Obligation. Proof. destruct AB as ((A & B) & HAB). simpl. t. Qed.
+
+    Program Definition tplu_co_z_points (P : point) (HPaff : z_of P = 1) : co_z_points :=
+      tplu P _.
+    Next Obligation. Proof. t. Qed.
+
+    Lemma opp_of_affine (P : Wpoint) (HPnz : P <> ∞ :> Wpoint) :
+      z_of (opp (of_affine P)) = 1.
+    Proof. t. Qed.
+
+    (* Scalar Multiplication on Weierstraß Elliptic Curves from Co-Z Arithmetic *)
+    (* Goundar, Joye, Miyaji, Rivain, Vanelli *)
+    (* Algorithm 14 Joye’s double-add algorithm with Co-Z addition formulæ *)
+    (* This is an adapted version that consumes and returns points in jacobian
+       coordinates, correctness assumes the scalar is odd (i.e., testbit 0 = true). *)
+    Definition joye_ladder_inner (scalarbitsz : Z) (testbit : Z -> bool) (P : point)
+      (HPaff : z_of P = 1) : point :=
+      (* Initialization *)
+      let b := testbit 1%Z in
+      let R1R0 := cswap_co_z_points b (tplu_co_z_points P HPaff) in
+      (* loop *)
+      let '(R1R0, _) :=
+        (@while (co_z_points*Z) (fun '(_, i) => (Z.ltb i scalarbitsz))
+           (fun '(R1R0, i) =>
+              let b := testbit i in
+              let R1R0 := cswap_co_z_points b R1R0 in
+              let R1R0 := cswap_co_z_points b (zdau_co_z_points R1R0) in
+              let i := Z.succ i in
+              (R1R0, i))
+           (Z.to_nat scalarbitsz) (* bound on loop iterations *)
+           (R1R0, 2%Z)) in
+      snd (proj1_sig R1R0).
+
+    (* Wrapper around joye_ladder_inner for points in affine coordinates,
+       it also readjusts the result if the scalar input is even. *)
+    Program Definition joye_ladder (scalarbitsz : Z) (testbit : Z -> bool) (P : Wpoint)
+      (HPnz : P <> ∞ :> Wpoint) : Wpoint :=
+      to_affine (
+      let P := of_affine P in
+      let R0 := joye_ladder_inner scalarbitsz testbit P _ in
+      let b := testbit 0%Z in
+      let mP := opp P in
+      (* Make sure R0 and -P are co-z *)
+      let R0R1 := make_co_z_points R0 mP (opp_of_affine _ HPnz) in
+      (* if b = 0 then R0 <- R0 - P and R1 <- R0 *)
+      (* if b = 1 then R0 <- R0 and R1 <- R0 - P *)
+      let R0 := fst (proj1_sig (cswap_co_z_points b (zaddu_co_z_points R0R1))) in
+      R0).
+    Next Obligation. Proof. t. Qed.
+
+    Section Proofs.
+
+      Local Notation scalarmult := (@scalarmult_ref Wpoint Wadd Wzero Wopp).
+      Local Notation scalarmult':= (@scalarmult_ref point add zero opp).
+
+      Local Instance Equivalence_Weq : Equivalence Weq.
+      Proof. apply Wgroup.(Hierarchy.commutative_group_group).(Hierarchy.group_monoid).(Hierarchy.monoid_Equivalence). Qed.
+
+      Lemma Pgroup :
+        @Hierarchy.group point eq add zero opp.
+      Proof.
+        destruct (@Group.group_by_isomorphism _ _ _ _ _ point eq add zero opp of_affine to_affine Wgroup.(Hierarchy.commutative_group_group) ltac:(apply Jacobian.to_affine_of_affine)); auto.
+        - intros; split; intros.
+          + rewrite <- (Jacobian.of_affine_to_affine a0), <- (Jacobian.of_affine_to_affine b0).
+            rewrite H. reflexivity.
+          + apply Jacobian.Proper_to_affine; auto.
+        - apply Jacobian.to_affine_add.
+        - intros. destruct a0 as (((X & Y) & Z) & HP).
+          unfold to_affine, Wopp, opp, Weq. simpl.
+          t.
+        - rewrite Jacobian.to_affine_of_affine. reflexivity.
+      Qed.
+
+      Lemma scalarmult_scalarmult' (n : Z) (P : Wpoint) :
+        eq (of_affine (scalarmult n P)) (scalarmult' n (of_affine P)).
+      Proof.
+        eapply (@homomorphism_scalarmult Wpoint Weq Wadd Wzero Wopp Wgroup.(Hierarchy.commutative_group_group) point eq add zero opp Pgroup scalarmult (@scalarmult_ref_is_scalarmult Wpoint Weq Wadd Wzero Wopp Wgroup.(Hierarchy.commutative_group_group)) scalarmult' (@scalarmult_ref_is_scalarmult point eq add zero opp Pgroup) of_affine ltac:(econstructor; [eapply Jacobian.of_affine_add|eapply Jacobian.Proper_of_affine])).
+      Qed.
+
+      Context {scalarbitsz : Z}
+              {Hscalarbitsz : (2 <= scalarbitsz)%Z}
+              {ordP : Z} {HordPpos : (2 < ordP)%Z}
+              {HordPodd : Z.odd ordP = true :> bool}.
+
+      Section Inner.
+        (* Proofs about joye_ladder_inner *)
+
+        (* Bit 0 of the scalar input is irrelevant *)
+        Lemma joye_ladder_inner_bit0_irr (bitsz : Z) (testbit0 testbit1 : Z -> bool)
+          (P : point) (HPaff : z_of P = 1)
+          (bit0_irr : forall i, (i >= 1)%Z -> testbit0 i = testbit1 i :> bool) :
+          eq (joye_ladder_inner bitsz testbit0 P HPaff)
+             (joye_ladder_inner bitsz testbit1 P HPaff).
+        Proof.
+          unfold joye_ladder_inner.
+          rewrite (surjective_pairing (while _ _ _ (cswap_co_z_points (testbit0 _) _, _))).
+          rewrite (surjective_pairing (while _ _ _ (cswap_co_z_points (testbit1 _) _, _))).
+          rewrite bit0_irr by lia.
+          match goal with
+          | |- eq (snd (proj1_sig (fst (while ?T0 ?B0 ?F0 ?I0))))
+                 (snd (proj1_sig (fst (while ?T1 ?B1 ?F1 ?I1)))) =>
+              set (test0 := T0);
+              set (body0 := B0);
+              set (fuel := F0);
+              set (init0 := I0);
+              set (body1 := B1)
+          end.
+          apply (while.preservation test0 body0 test0 body1 (fun s1 s2 => eq (snd (proj1_sig (fst s1))) (snd (proj1_sig (fst s2))) /\ (s1 = s2 :> (co_z_points * Z)) /\ (2 <= snd s1)%Z)).
+          - intros s1 s2 (_ & <- & _). reflexivity.
+          - unfold test0. intros (PQ1 & i1) (PQ2 & i2) _.
+            cbn [fst snd]. intros (_ & Heq & Hi). inversion Heq; clear Heq.
+            subst PQ1; subst i1.
+            unfold body0, body1. rewrite bit0_irr by lia.
+            split; [reflexivity|split; [reflexivity|simpl; lia] ].
+          - repeat (split; try reflexivity).
+        Qed.
+
+        Context {n : Z} {Hnodd : n = Z.setbit n 0 :> Z}
+          {Hn : (2 <= n < 2^scalarbitsz)%Z}
+          {P : point} {HPaff : z_of P = 1}
+          {HordP : forall l, (eq (scalarmult' l P) zero) <-> exists k, (l = k * ordP :> Z)%Z}.
+        Local Notation testbitn := (Z.testbit n).
+        Context {HSS : forall i, (2 <= i <= scalarbitsz)%Z -> not (eq (scalarmult' (SS n (Z.to_nat i)) P) zero)}
+                {HTT : forall i, (2 <= i <= scalarbitsz)%Z -> not (eq (scalarmult' (TT n (Z.to_nat i)) P) zero)}.
+
+        Lemma mult_two_power (k : Z) :
+          (0 <= k)%Z ->
+          not (eq (scalarmult' (2^k)%Z P) zero).
+        Proof.
+          eapply (Zlt_0_ind (fun x => ~ eq (scalarmult' (2 ^ x) P) zero)).
+          intros x Hind Hxpos Heqz.
+          destruct (proj1 (HordP (2^x)%Z) Heqz) as [l Hl].
+          destruct (Z.eq_dec x 0); [subst x|].
+          - simpl in Hl. clear -Hl HordPpos.
+            generalize (Z.divide_1_r_nonneg ordP ltac:(lia) ltac:(exists l; lia)).
+            lia.
+          - assert (x = Z.succ (Z.pred x) :> Z) by lia.
+            rewrite H in Hl. rewrite Z.pow_succ_r in Hl; [|lia].
+            generalize (Znumtheory.prime_mult 2%Z Znumtheory.prime_2 l ordP ltac:(exists (2 ^ Z.pred x)%Z; lia)).
+            intros [A|A]; destruct A as [m Hm]; [|replace ordP with (0 + 2 * m)%Z in HordPodd by lia; rewrite Z.odd_add_mul_2 in HordPodd; simpl in HordPodd; congruence].
+            rewrite Hm in Hl.
+            assert ((2 ^ Z.pred x)%Z = (m * ordP)%Z :> Z) by lia.
+            apply (Hind (Z.pred x) ltac:(lia)).
+            eapply HordP. exists m; assumption.
+        Qed.
+
+        Lemma mult_two (k : Z) :
+          eq (scalarmult' (2 * k)%Z P) zero ->
+          eq (scalarmult' k P) zero.
+        Proof.
+          intros Heqz. destruct (proj1 (HordP (2 * k)%Z) Heqz) as [l Hl].
+          generalize (Znumtheory.prime_mult 2%Z Znumtheory.prime_2 l ordP ltac:(exists k; lia)).
+          intros [A|A]; destruct A as [m Hm]; [|replace ordP with (0 + 2 * m)%Z in HordPodd by lia; rewrite Z.odd_add_mul_2 in HordPodd; simpl in HordPodd; congruence].
+          rewrite Hm in Hl. assert (k = m * ordP :> Z)%Z as -> by lia.
+          apply HordP; eauto.
+        Qed.
+
+        Lemma HSS_plus_TT (m : Z) (k : nat) :
+          not (eq (scalarmult' (SS m k + TT m k)%Z P) zero).
+        Proof. rewrite SS_plus_TT. apply mult_two_power. lia. Qed.
+
+        Lemma SS1 : SS n 1 = 1%Z :> Z.
+        Proof. cbv [SS]. rewrite Z.pow_1_r, <-Z.bit0_mod, Hnodd, Z.setbit_eqb; trivial; lia. Qed.
+
+        Lemma TT1 : TT n 1 = 1%Z :> Z.
+        Proof. cbv [TT]. rewrite SS1; trivial. Qed.
+
+        Lemma SS2 : SS n 2 = (if Z.testbit n 1 then 3%Z else 1%Z) :> Z.
+        Proof.
+          cbv [SS]; rewrite Hnodd, <-Z.land_ones, Z.land_comm by lia.
+          destruct n; cbn; trivial; repeat (destruct p; cbn; trivial).
+        Qed.
+
+        Lemma TT2 : TT n 2 = (if Z.testbit n 1 then 1%Z else 3%Z) :> Z.
+        Proof. cbv [TT]. rewrite SS2. case Z.testbit; lia. Qed.
+
+        Lemma SS_TT2 : (if testbitn 1 then SS n 2 else TT n 2) = 3%Z :> Z.
+        Proof. rewrite SS2, TT2; try lia; case testbitn; trivial. Qed.
+
+        Lemma HordP3 :
+          not (eq (scalarmult' 3%Z P) zero).
+        Proof.
+          rewrite <- SS_TT2, <- (Nat2Z.id 2).
+          case testbitn; [eapply HSS|eapply HTT]; lia.
+        Qed.
+
+        Hint Unfold fst snd proj1_sig : points_as_coordinates.
+        Hint Unfold fieldwise fieldwise' : points_as_coordinates.
+
+        Lemma eq_fieldwise (P1 P2 : point) :
+          fieldwise (n:=3) Feq
+            (proj1_sig P1) (proj1_sig P2) ->
+          eq P1 P2.
+        Proof. clear -field; intros; t. Qed.
+
+        Lemma Pynz :
+          y_of P <> 0.
+        Proof.
+          intro Hy. assert (HA : eq P (opp P)).
+          - apply eq_fieldwise. destruct P as (((X & Y) & Z) & HP).
+            simpl; cbv in HPaff, Hy; repeat split; fsatz.
+          - apply (mult_two_power 1%Z ltac:(lia)).
+            replace (2 ^ 1)%Z with (1 - -1)%Z by lia.
+            rewrite (scalarmult_sub_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+            rewrite <- (scalarmult_1_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))) in HA.
+            rewrite HA.
+            rewrite <- (scalarmult_opp_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+            rewrite <- (scalarmult_sub_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+            replace (- (1) - -1)%Z with 0%Z by lia. reflexivity.
+        Qed.
+
+        Lemma add_opp_zero (A : point) :
+          eq (add A (opp A)) zero.
+        Proof.
+          generalize (Jacobian.add_opp A).
+          destruct (add A (opp A)) as (((X & Y) & Z) & H).
+          cbn. intros HP; destruct (dec (Z = 0)); fsatz.
+        Qed.
+
+        Lemma scalarmult_difference (A : point) (n1 n2 : Z):
+          eq (scalarmult' n1 A) (scalarmult' n2 A) ->
+          eq (scalarmult' (n1 - n2)%Z A) zero.
+        Proof.
+          intros. rewrite (scalarmult_sub_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))), H, <- (scalarmult_sub_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))), Z.sub_diag.
+          apply (scalarmult_0_l (is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+        Qed.
+
+        Lemma dblu_scalarmult' (Q : point) (Hz1 : z_of Q = 1) (Hynz : y_of Q <> 0) :
+          let '(M, N) := dblu Q Hz1 in
+          eq M (scalarmult' 2 Q)
+          /\ eq N (scalarmult' 1 Q).
+        Proof.
+          generalize (@Jacobian.dblu_correct F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field char_ge_3 Feq_dec char_ge_12 Q Hz1 Hynz).
+          rewrite (surjective_pairing (dblu _ _)) at 1.
+          intros (A & B & C). split.
+          - rewrite <- A. replace 2%Z with (1 + 1)%Z by lia.
+            rewrite (scalarmult_add_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+            rewrite (scalarmult_1_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+            rewrite <- Jacobian.add_double; reflexivity.
+          - rewrite (scalarmult_1_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+            symmetry. assumption.
+        Qed.
+
+        Lemma co_xz_implies (P1 P2 : point) (Hxeq : x_of P1 = x_of P2)
+          (Hzeq : z_of P1 = z_of P2) :
+          eq P1 P2 \/ eq P1 (opp P2).
+        Proof.
+          clear -Hxeq Hzeq. prept; [tauto|fsatz|fsatz|].
+          assert (f4 = f1 \/ f4 = Fopp f1) by (destruct (dec (f4 = f1)); [left; assumption|right; fsatz]).
+          destruct H; [left; repeat split; fsatz|right; repeat split; fsatz].
+        Qed.
+
+        Lemma tplu_scalarmult' {p q} (H : tplu P HPaff = (p, q) :> _) :
+          eq p (scalarmult' 3 P) /\ eq q (scalarmult' 1 P) /\ co_z p q.
+        Proof.
+          intros; unfold tplu. generalize (dblu_scalarmult' P HPaff Pynz).
+          inversion_prod; subst p q.
+          rewrite (surjective_pairing (dblu _ _)) at 1.
+          set (P1 := (snd (dblu P HPaff))).
+          set (P2 := (fst (dblu P HPaff))). intros [A1 A2].
+          destruct (dec (x_of P1 = x_of P2)) as [Hxe|Hxne].
+          { destruct (co_xz_implies P1 P2 Hxe (CoZ.Jacobian.tplu_obligation_1 P HPaff)) as [A|A].
+            - rewrite A1, A2 in A. elim (mult_two_power 1%Z ltac:(lia)).
+              rewrite <- A. replace 1%Z with (2 - 1)%Z by lia.
+              apply scalarmult_difference; symmetry; assumption.
+            - rewrite A1, A2, <- (scalarmult_opp_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))) in A.
+              apply scalarmult_difference in A.
+              elim HordP3. exact A. }
+          generalize (Jacobian.zaddu_correct _ _ (CoZ.Jacobian.tplu_obligation_1 P HPaff) Hxne).
+          rewrite (surjective_pairing (zaddu _ _ _)) at 1.
+          intros (A & B & C); subst P1 P2.
+          repeat try split; trivial.
+          { rewrite <-A, A1, A2, (@scalarmult_add_l point eq add zero opp Pgroup scalarmult' (@scalarmult_ref_is_scalarmult _ _ _ _ _ Pgroup) 1 2); reflexivity. }
+          { rewrite <-B. rewrite A2. reflexivity. }
+        Qed.
+
+        (* Since co-Z formulas are incomplete, we need to show that we won't hit the neutral point for ZDAU in the loop *)
+        Lemma zaddu_SS_TT (i : Z) (B1 B2 Y1 Y2 R0 R1 : point) (HB12 : co_z B1 B2)
+          (Hi : (2 <= i < scalarbitsz)%Z)
+          (HBx : x_of B1 <> x_of B2)
+          (HY : zaddu B1 B2 HB12 = (Y1, Y2) :> point * point)
+          (HR0 : eq R0 (scalarmult' (SS n (Z.to_nat i)) P))
+          (HR1 : eq R1 (scalarmult' (TT n (Z.to_nat i)) P))
+          (HB1 : B1 = (if testbitn i then R0 else R1) :> point)
+          (HB2 : B2 = (if testbitn i then R1 else R0) :> point) :
+          x_of Y1 <> x_of Y2.
+        Proof.
+          intro XX. generalize (Jacobian.zaddu_correct B1 B2 HB12 HBx).
+          rewrite HY. intros (W1 & W2 & W3).
+          destruct (co_xz_implies _ _ XX W3) as [W|W]; rewrite W, <- W2, HB1, HB2 in W1.
+          - destruct (testbitn i);
+            rewrite HR1, HR0, <- (scalarmult_add_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))) in W1;
+            apply scalarmult_difference in W1;
+            rewrite Z.add_simpl_l in W1;
+            [apply (HTT i ltac:(lia))|apply (HSS i ltac:(lia))]; auto.
+          - destruct (testbitn i) eqn:Hti;
+            rewrite HR1, HR0, <- (scalarmult_add_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))), <- (scalarmult_opp_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))) in W1;
+            apply scalarmult_difference in W1.
+            all: match goal with
+                 | H : eq (scalarmult' ?X _) zero |- _ =>
+                     match X with
+                     | (SS _ _ + _ - _)%Z =>
+                         replace X with (SS n (S (Z.to_nat i))) in H by (rewrite SS_succ, Z2Nat.id, Hti; lia)
+                     | (TT _ _ + _ - _)%Z =>
+                         replace X with (TT n (S (Z.to_nat i))) in H by (rewrite TT_succ, Z2Nat.id, Hti; lia)
+                     end
+                 end.
+            all: rewrite <- Z2Nat.inj_succ in W1; try lia.
+            all: match goal with
+                 | H : eq (scalarmult' (SS _ _) _) zero |- _ =>
+                     apply (HSS (Z.succ i) ltac:(lia))
+                 | _ =>
+                     apply (HTT (Z.succ i) ltac:(lia))
+                 end; auto.
+        Qed.
+
+        Lemma SS_TT_xne (i : Z) (R0 R1 : point) (HCOZ : co_z R0 R1)
+          (Hi : (2 <= i < scalarbitsz)%Z)
+          (HR0 : eq R0 (scalarmult' (SS n (Z.to_nat i)) P))
+          (HR1 : eq R1 (scalarmult' (TT n (Z.to_nat i)) P)) :
+          x_of R0 <> x_of R1.
+        Proof.
+          assert (HH : forall A, eq (opp A) zero -> eq A zero) by (clear -field; unfold zero, Wzero; intros; t).
+          intros Hxe. generalize (co_xz_implies _ _ Hxe HCOZ).
+          destruct (Z.eq_dec 2 i).
+          { subst i. replace (Z.to_nat 2) with 2%nat in HR1 by lia.
+            replace (Z.to_nat 2) with 2%nat in HR0 by lia.
+            rewrite TT2 in HR1. rewrite SS2 in HR0.
+            rewrite HR0, HR1, <- (scalarmult_opp_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+            destruct (testbitn 1) eqn:Htj; intros [Q|Q]; apply scalarmult_difference in Q.
+            all: match goal with
+                 | H : eq (scalarmult' ?X _) zero |- _ =>
+                     try replace X with 2%Z in H by lia;
+                     try replace X with 4%Z in H by lia;
+                     try replace X with (- (2))%Z in H by lia
+                 end.
+            all: match goal with
+                 | H : eq (scalarmult' ?X _) zero |- _ =>
+                     match X with
+                     | 2%Z => apply (mult_two_power 1%Z ltac:(lia) H)
+                     | 4%Z => apply (mult_two_power 2%Z ltac:(lia) H)
+                     | (- (2))%Z =>
+                         rewrite (scalarmult_opp_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))) in Q;
+                         apply HH in H;
+                         apply (mult_two_power 1%Z ltac:(lia) H)
+                     end
+                 end. }
+          { set (j := Z.pred i). assert (His : i = Z.succ j :> Z) by lia.
+            assert (Hj : (2 <= j)%Z) by lia.
+            rewrite His, Z2Nat.inj_succ in HR0, HR1 by lia.
+            rewrite HR0, HR1, <- (scalarmult_opp_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+            rewrite TT_succ, SS_succ, Z2Nat.id by lia.
+            destruct (testbitn j) eqn:Htj; intros [Q|Q]; apply scalarmult_difference in Q.
+            all: try match goal with
+                   | H : eq (scalarmult' ?X _) zero |- _ =>
+                       match X with
+                       | Z.sub _ ?Y =>
+                           match Y with
+                           | (- _)%Z =>
+                               replace X with (2 * (SS n (Z.to_nat j) + TT n (Z.to_nat j)))%Z in H by lia;
+                               apply mult_two in H;
+                               apply (HSS_plus_TT n (Z.to_nat j) H)
+                           | (TT _ _) =>
+                               replace X with (2 * SS n (Z.to_nat j))%Z in H by lia;
+                               apply mult_two in H;
+                               apply (HSS j ltac:(lia) H)
+                           | (SS _ _) =>
+                               replace X with (2 * TT n (Z.to_nat j))%Z in H by lia;
+                               apply mult_two in H;
+                               apply (HTT j ltac:(lia) H)
+                           | (Z.add _ _) =>
+                               replace X with (2 * - TT n (Z.to_nat j))%Z in H by lia;
+                               apply mult_two in H;
+                               rewrite (scalarmult_opp_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))) in H;
+                               apply HH in H;
+                               apply (HTT j ltac:(lia) H)
+                           end
+                       end
+                   end. }
+        Qed.
+
+        (* When n is odd, joye_ladder_inner computes [n]P *)
+        Lemma joye_ladder_inner_correct :
+          eq (joye_ladder_inner scalarbitsz testbitn P HPaff)
+             (scalarmult' n P).
+        Proof.
+          unfold joye_ladder_inner. set (WW := while _ _ _ _).
+          destruct (tplu_co_z_points P HPaff) as ((A1 & A2) & HA12) eqn:Htplu.
+          rewrite (sig_eta (tplu_co_z_points _ _)) in Htplu.
+          apply proj1_sig_eq in Htplu; simpl in Htplu.
+          (* Initialize the ladder state with ([3]P, [1]P) or its symmetric *)
+          destruct (cswap_co_z_points (testbitn 1) _) as ((A3 & A4) & HA34) eqn:HA1.
+          rewrite (sig_eta (cswap_co_z_points _ _)) in HA1.
+          apply proj1_sig_eq in HA1. cbn [proj1_sig cswap_co_z_points] in HA1.
+          assert (A3 = (if testbitn 1 then A2 else A1) :> point) as HA3 by (destruct (testbitn 1); inversion HA1; auto).
+          assert (A4 = (if testbitn 1 then A1 else A2) :> point) as HA4 by (destruct (testbitn 1); inversion HA1; auto).
+          clear HA1. destruct (tplu_scalarmult' Htplu) as (HeqA1 & HeqA2 & _).
+          set (inv :=
+                 fun '(R1R0, i) =>
+                   let '(R1, R0) := proj1_sig (R1R0:co_z_points) in
+                   (2 <= i <= scalarbitsz)%Z /\
+                     (eq R1 (scalarmult' (TT n (Z.to_nat i)) P)
+                      /\ eq R0 (scalarmult' (SS n (Z.to_nat i)) P))
+                   /\ ((i < scalarbitsz)%Z -> x_of R1 <> x_of R0)).
+          assert (HH : forall (A B : Prop), A -> (A -> B) -> A /\ B) by tauto.
+          assert (WWinv : inv WW /\ (snd WW = scalarbitsz :> Z)).
+          { set (measure := fun (state : (co_z_points*Z)) => ((Z.to_nat scalarbitsz) + 2 - (Z.to_nat (snd state)))%nat).
+            unfold WW. replace (Z.to_nat scalarbitsz) with (measure (exist _ (A3, A4) HA34, 2%Z)) by (unfold measure; simpl; lia).
+            eapply (while.by_invariant inv measure (fun s => inv s /\ (snd s = scalarbitsz :> Z))).
+            - (* Invariant holds at beginning *)
+              unfold inv. cbn [proj1_sig].
+              split; [lia|]. apply HH.
+              + change (Z.to_nat 2) with 2%nat.
+                rewrite SS2, TT2, HA3, HA4.
+                case Z.testbit; auto.
+              + intros [He1 He2] Hxe. symmetry.
+                apply (SS_TT_xne 2%Z A4 A3 ltac:(apply co_z_comm; exact HA34) ltac:(lia)); eauto.
+            - (* Invariant is preserved by the loop,
+               measure decreases,
+               and post-condition i = scalarbitsz. *)
+              intros s Hs. destruct s as (R1R0 & i).
+              destruct R1R0 as ((R1 & R0) & HCOZ).
+              destruct Hs as (Hi & (HR1 & HR0) & Hx).
+              destruct (Z.ltb i scalarbitsz) eqn:Hltb.
+              + apply Z.ltb_lt in Hltb.
+                split.
+                * (* Invariant preservation.
+                   The loop body is basically :
+                   (R1, R0) <- cswap (testbitn i) (R1, R0);
+                   (R1, R0) <- ZDAU (R1, R0);
+                   (R1, R0) <- cswap (testbitn i) (R1, R0);
+                   *)
+                  (* Start by giving names to all intermediate values *)
+                  unfold inv. destruct (cswap_co_z_points (testbitn i) (exist _ _ _)) as ((B1 & B2) & HB12) eqn:Hswap1.
+                  rewrite (sig_eta (cswap_co_z_points _ _)) in Hswap1.
+                  apply proj1_sig_eq in Hswap1. simpl in Hswap1.
+                  assert (HB1: B1 = (if testbitn i then R0 else R1) :> point) by (destruct (testbitn i); congruence).
+                  assert (HB2: B2 = (if testbitn i then R1 else R0) :> point) by (destruct (testbitn i); congruence).
+                  clear Hswap1.
+                  destruct (zdau_co_z_points _) as ((C1 & C2) & HC12) eqn:HZDAU.
+                  rewrite (sig_eta (zdau_co_z_points _)) in HZDAU.
+                  apply proj1_sig_eq in HZDAU. simpl in HZDAU.
+                  assert (HBx : x_of B1 <> x_of B2) by (rewrite HB1, HB2; destruct (testbitn i); [symmetry|]; auto).
+                  destruct (zaddu B1 B2 (zdau_co_z_points_obligation_1 (exist (fun '(A, B) => co_z A B) (B1, B2) HB12) B1 B2 eq_refl)) as (Y1 & Y2) eqn:HY.
+                  assert (HYx : x_of Y1 <> x_of Y2) by (eapply zaddu_SS_TT; eauto; lia).
+                  generalize (@Jacobian.zdau_correct_alt F Feq Fzero Fone Fopp Fadd Fsub Fmul Finv Fdiv a b field char_ge_3 Feq_dec char_ge_12 ltac:(unfold id in *; fsatz) B1 B2 (zdau_co_z_points_obligation_1 (exist (fun '(A, B) => co_z A B) (B1, B2) HB12) B1 B2 eq_refl) HBx ltac:(rewrite HY; simpl; apply HYx)).
+                  rewrite HZDAU. intros (HC1 & HC2 & _).
+                  destruct (cswap_co_z_points (testbitn i) _) as ((D1 & D2) & HD12) eqn:HD.
+                  rewrite (sig_eta (cswap_co_z_points _ _)) in HD.
+                  apply proj1_sig_eq in HD. cbn [proj1_sig cswap_co_z_points] in HD.
+                  assert (HD1 : D1 = (if testbitn i then C2 else C1) :> point) by (destruct (testbitn i); congruence).
+                  assert (HD2 : D2 = (if testbitn i then C1 else C2) :> point) by (destruct (testbitn i); congruence).
+                  clear HD. simpl.
+                  (* invariant preservation *)
+                  (* counter still within bounds *)
+                  split; [lia|]. rewrite HD1, HD2. apply HH.
+                  { (* New values are indeed [SS (i+1)]P and [TT (i+1)]P *)
+                    destruct (testbitn i) eqn:Hti;
+                    rewrite <- HC1, <- HC2, HB1, HB2;
+                    replace (Z.to_nat (Z.succ i)) with (S (Z.to_nat i)) by lia;
+                    rewrite SS_succ, TT_succ, Z2Nat.id by lia;
+                    rewrite Hti; split; try assumption;
+                    rewrite <- Jacobian.add_double; try reflexivity;
+                    rewrite HR0, HR1;
+                    repeat rewrite <- (scalarmult_add_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup)));
+                    rewrite <- Z.add_diag; reflexivity. }
+                  { (* Make sure we don't hit bad values *)
+                    intros [He1 He2] Hxe. symmetry. eapply (SS_TT_xne (Z.succ i)); eauto.
+                    - destruct (testbitn i); [|apply co_z_comm]; auto.
+                    - lia. }
+                * (* measure decreases *)
+                  apply Z.ltb_lt in Hltb.
+                  unfold measure; simpl; lia.
+              + (* Post-condition *)
+                simpl; split; auto.
+                rewrite Z.ltb_nlt in Hltb. lia. }
+          destruct WWinv as (Hinv & Hj).
+          destruct WW as (R1R0 & i). destruct R1R0 as ((R1 & R0) & HCOZ).
+          simpl in Hj; subst i. destruct Hinv as (_ & (_ & HR0) & _).
+          rewrite (SSn n scalarbitsz ltac:(lia) ltac:(lia)) in HR0.
+          exact HR0.
+        Qed.
+      End Inner.
+
+      (* We compute [n]P where P ≠ ∞ and n < ord(P) *)
+      Context {n : Z} {Hn : (2 <= n < 2^scalarbitsz)%Z}
+              {P : Wpoint} {HPnz : P <> ∞ :> Wpoint}
+              {HordP : forall l, (Weq (scalarmult l P) ∞) <-> exists k, (l = k * ordP :> Z)%Z}
+              {HordPn : (n + 2 < ordP)%Z}.
+
+      Local Notation testbitn := (Z.testbit n).
+      (* Local Notation n' := (if testbitn 0 then n else (n + 1)%Z). *)
+      Local Notation n' := (Z.setbit n 0) (only parsing).
+      Local Notation testbitn' := (Z.testbit n') (only parsing).
+
+      (* Technically, if q is the order of the curve, and scalarbitsz = [log2 q],
+         and ordP is close to q, then only the last two values of SS and TT are
+         dangerous.
+         See §3.1 of "Faster Montgomery and double-add ladders for short
+         Weierstrass curves" (Mike Hamburg, CHES'20) for instance.
+       *)
+      Context {HSS : forall i, (2 <= i <= scalarbitsz)%Z -> not (Weq (scalarmult (SS n' (Z.to_nat i)) P) ∞) }
+              {HTT : forall i, (2 <= i <= scalarbitsz)%Z -> not (Weq (scalarmult (TT n' (Z.to_nat i)) P) ∞) }.
+
+      Lemma Hn' : (2 <= n' < 2^scalarbitsz)%Z.
+      Proof.
+        rewrite Z.setbit_spec'; simpl; split.
+        - etransitivity; [|apply LandLorShiftBounds.Z.lor_lower]; lia.
+        - apply (LandLorShiftBounds.Z.lor_range n 1 scalarbitsz); lia.
+      Qed.
+
+      Lemma scalarmult_eq_weq_conversion (k1 k2 : Z) :
+        Weq (scalarmult k1 P) (scalarmult k2 P) <-> eq (scalarmult' k1 (of_affine P)) (scalarmult' k2 (of_affine P)).
+      Proof.
+        split; intros.
+        - repeat rewrite <- scalarmult_scalarmult'.
+          eapply Jacobian.Proper_of_affine. apply H.
+        - rewrite <- Jacobian.to_affine_of_affine at 1.
+          symmetry. rewrite <- Jacobian.to_affine_of_affine at 1.
+          apply Jacobian.Proper_to_affine.
+          symmetry; repeat rewrite scalarmult_scalarmult'; auto.
+      Qed.
+
+      Lemma HordP_alt (k : Z) :
+        (0 < k < ordP)%Z ->
+        not (Weq (scalarmult k P) ∞).
+      Proof.
+        intros Hbds Heq. destruct (proj1 (HordP k) Heq) as [l Hl].
+        clear -Hbds Hl. generalize (Zmult_gt_0_lt_0_reg_r ordP l ltac:(lia) ltac:(lia)).
+        intros. generalize (proj1 (Z.mul_le_mono_pos_r 1%Z l ordP ltac:(lia)) ltac:(lia)). lia.
+      Qed.
+
+      Lemma n'_alt :
+        n' = (if testbitn 0 then n else (n + 1)%Z) :> Z.
+      Proof.
+        apply Z.bits_inj'; intros.
+        destruct (Z.eq_dec n0 0) as [->|?]; [rewrite Z.setbit_eq|rewrite Z.setbit_neq]; try lia.
+        - destruct (testbitn 0) eqn:Hbit0; auto.
+          rewrite Z.bit0_odd, <- Z.even_pred.
+          replace (Z.pred (n + 1))%Z with n by lia.
+          rewrite <- Z.negb_odd, <- Z.bit0_odd, Hbit0; reflexivity.
+        - destruct (testbitn 0) eqn:Hbit0; auto.
+          replace n0 with (Z.succ (n0 - 1))%Z by lia.
+          rewrite Z.bit0_odd in Hbit0.
+          rewrite (Zeven_div2 n); [|apply Zeven_bool_iff; rewrite <- Z.negb_odd, Hbit0; reflexivity].
+          rewrite Z.testbit_even_succ, Z.testbit_odd_succ; auto; lia.
+      Qed.
+
+      Lemma joye_ladder_correct :
+        Weq (joye_ladder scalarbitsz testbitn P HPnz) (scalarmult n P).
+      Proof.
+        (* Unfold the ladder *)
+        rewrite <- (Jacobian.to_affine_of_affine (scalarmult n P)).
+        apply Jacobian.Proper_to_affine. rewrite scalarmult_scalarmult'.
+        set (P' := of_affine P). set (HPaff := joye_ladder_obligation_1 P HPnz).
+        assert (Hnodd : n' = Z.setbit n' 0 :> Z) by (repeat rewrite Z.setbit_spec'; rewrite <- Z.lor_assoc, Z.lor_diag; reflexivity).
+        assert (HordP' : forall l, (eq (scalarmult' l P') zero) <-> exists k, (l = k * ordP :> Z)%Z).
+        { intros l; split; replace zero with (scalarmult' 0%Z P') by reflexivity.
+          - intros H; apply scalarmult_eq_weq_conversion in H; apply HordP; auto.
+          - intros H; apply scalarmult_eq_weq_conversion; apply HordP; auto. }
+        assert (HSS' : forall i, (2 <= i <= scalarbitsz)%Z -> not (eq (scalarmult' (SS n' (Z.to_nat i)) P') zero)).
+        { replace zero with (scalarmult' 0%Z P') by reflexivity.
+          intros i Hi He; apply scalarmult_eq_weq_conversion in He; apply (HSS i Hi); auto. }
+        assert (HTT' : forall i, (2 <= i <= scalarbitsz)%Z -> not (eq (scalarmult' (TT n' (Z.to_nat i)) P') zero)).
+        { replace zero with (scalarmult' 0%Z P') by reflexivity.
+          intros i Hi He; apply scalarmult_eq_weq_conversion in He; apply (HTT i Hi); auto. }
+        generalize (joye_ladder_inner_correct (n:=n') (P:=of_affine P) (HPaff:=HPaff) (Hnodd:=Hnodd) (Hn:=Hn') (HordP:=HordP')(HSS:=HSS') (HTT:=HTT')).
+        intros Hinner. rewrite (joye_ladder_inner_bit0_irr scalarbitsz testbitn' testbitn (of_affine P) HPaff ltac:(intros; rewrite Z.setbit_neq; trivial; lia)) in Hinner.
+        set (R0 := joye_ladder_inner scalarbitsz testbitn P' HPaff).
+        cbv zeta. fold R0. fold P' in Hinner. fold R0 in Hinner.
+        (* We now have R0 = [n']P *)
+        destruct (make_co_z_points R0 _ _) as ((B1 & B2) & HB12) eqn:HMCZ.
+        rewrite (sig_eta (make_co_z_points _ _ _)) in HMCZ.
+        apply proj1_sig_eq in HMCZ; simpl in HMCZ.
+        (* Prove [n']P ≠ ∞ *)
+        assert (HR0znz : z_of R0 <> 0).
+        { intro. apply (HSS scalarbitsz ltac:(lia)).
+          replace Wzero with (scalarmult 0%Z P) by reflexivity.
+          apply scalarmult_eq_weq_conversion.
+          rewrite SSn by (generalize Hn'; lia).
+          fold P'. rewrite <- Hinner. destruct R0 as (((? & ?) & ?) & ?).
+          cbn in H; cbn. clear -field H.
+          destruct (dec (f1 = 0)); fsatz. }
+        (* Have co-Z representations of [n']P and [-1]P *)
+        generalize (Jacobian.make_co_z_correct R0 (opp P') (opp_of_affine P HPnz) HR0znz).
+        rewrite HMCZ. intros (HB1 & HB2 & _).
+        destruct (zaddu_co_z_points _) as ((C1 & C2) & HC12) eqn:HZADDU.
+        rewrite (sig_eta (zaddu_co_z_points _)) in HZADDU.
+        apply proj1_sig_eq in HZADDU. simpl in HZADDU.
+        (* Ensure ZADDU doesn't hit the neutral point *)
+        assert (Hxne: x_of B1 <> x_of B2).
+        { intro Hxe. destruct (co_xz_implies _ _ Hxe HB12) as [HEq|HNeq]; [rewrite <- HB1, Hinner, <- HB2 in HEq|rewrite <- HB1, Hinner, <- HB2 in HNeq].
+          - rewrite <- (scalarmult_opp1_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))) in HEq.
+            apply scalarmult_difference in HEq.
+            apply (HordP_alt (n' - -1)%Z).
+            + rewrite n'_alt; destruct (testbitn 0); lia.
+            + replace Wzero with (scalarmult 0 P) by reflexivity.
+              apply scalarmult_eq_weq_conversion. auto.
+          - rewrite (Group.inv_inv (H:=Pgroup)) in HNeq.
+            rewrite <- (scalarmult_1_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup)) P') in HNeq at 2.
+            apply scalarmult_difference in HNeq.
+            apply (HordP_alt (n' - 1)%Z).
+            + rewrite n'_alt; destruct (testbitn 0); lia.
+            + replace Wzero with (scalarmult 0 P) by reflexivity.
+              apply scalarmult_eq_weq_conversion. auto. }
+        generalize (Jacobian.zaddu_correct B1 B2 (zaddu_co_z_points_obligation_1 (exist (fun '(A, B) => co_z A B) (B1, B2) HB12) B1 B2 eq_refl) Hxne).
+        rewrite HZADDU. intros (HC1 & HC2 & _).
+        rewrite <- HB1, <- HB2, Hinner in HC1. rewrite <- HB1, Hinner in HC2.
+        destruct (cswap_co_z_points (testbitn 0) _) as ((D1 & D2) & HD12) eqn:Hswap.
+        rewrite (sig_eta (cswap_co_z_points _ _)) in Hswap.
+        apply proj1_sig_eq in Hswap. cbn [proj1_sig cswap_co_z_points] in Hswap.
+        simpl. assert (D1 = (if testbitn 0 then C2 else C1) :> point) as -> by (destruct (testbitn 0); congruence).
+        rewrite n'_alt in HC1, HC2.
+        destruct (testbitn 0); [rewrite <- HC2; reflexivity|rewrite <- HC1].
+        rewrite <- (scalarmult_opp1_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))), <- (scalarmult_add_l (groupG:=Pgroup) (mul_is_scalarmult:=scalarmult_ref_is_scalarmult (groupG:=Pgroup))).
+        replace (n + 1 + -1)%Z with n by lia. reflexivity.
+      Qed.
+    End Proofs.
+  End JoyeLadder.
+End ScalarMult.