\section{Proving that X25519 in Coq matches the mathematical model}
\label{sec3-maths}

In this section we first present the work of Bartzia and Strub \cite{DBLP:conf/itp/BartziaS14} (\ref{Weierstrass}).
We extend it to support Montgomery curves (\ref{montgomery}) with homogeneous coordinates (\ref{projective}) and prove the correctness of the ladder (\ref{ladder}).

We then prove the Montgomery ladder computes
the x-coordinate of scalar multiplication over $\F{p^2}$
(Theorem 2.1 by Bernstein \cite{Ber06}) where $p$ is the prime $\p$.

\subsection{Formalization of Elliptic Curves}

We consider elliptic curves over a field $\K$. We assume that the
characteristic of $\K$ is neither 2 or 3.

\begin{definition}
Let a field $\K$, using an appropriate choice of coordinates, an elliptic curve $E$
is a plane cubic algebraic curve $E(x,y)$ defined by an equation of the form:
$$E : y^2 + a_1 xy + a_3 y = x^3 + a_2 x^2 + a_4 x + a_6$$
where the $a_i$'s are in \K\ and the curve has no singular point (\ie no cusps
or self-intersections). The set of points, written $E(\K)$, is formed by the
solutions $(x,y)$ of $E$ augmented by a distinguished point $\Oinf$ (called point at infinity):
$$E(\K) = \{(x,y) \in \K \times \K | E(x,y)\} \cup \{\Oinf\}$$
\end{definition}

\subsubsection{Weierstra{\ss} Curves}
\label{Weierstrass}
This equation $E(x,y)$ can be reduced into its short Weierstra{\ss} form.

\begin{definition}
Let $a \in \K$, and $b \in \K$ such that $$\Delta(a,b) = -16(4a^3 + 27b^2) \neq 0.$$
The \textit{elliptic curve} $E_{a,b}(\K)$ is the set of all points $(x,y) \in \K^2$ satisfying the equation:
$$y^2 = x^3 + ax + b,$$
along with an additional formal point $\Oinf$, ``at infinity''. Such curve does not present any singularity.
\end{definition}

In this setting, Bartzia and Strub defined the parametric type \texttt{ec} which
represent the points on a specific curve. It is parameterized by
a \texttt{K : ecuFieldType} -- the type of fields which characteristic is not 2 or 3 --
and \texttt{E : ecuType} -- a record that packs the curve parameters $a$ and $b$
along with the proof that $\Delta(a,b) \neq 0$.
\begin{lstlisting}[language=Coq]
Inductive point := EC_Inf | EC_In of K * K.
Notation "(| x, y |)" := (EC_In x y).
Notation "\infty" := (EC_Inf).

Record ecuType :=
  { A : K; B : K; _ : 4 * A^3 + 27 * B^2 != 0}.
Definition oncurve (p : point) :=
  if p is (| x, y |)
    then y^2 == x^3 + A * x + B
    else true.
Inductive ec : Type := EC p of oncurve p.
\end{lstlisting}

Points of an elliptic curve can be equipped with a structure of an abelian group.
\begin{itemize}
  \item The negation of a point $P = (x,y)$ by taking the symmetric with respect to the x axis $-P = (x, -y)$.
  \item The addition of two points $P$ and $Q$ is defined by the negation of third intersection
  of the line passing by $P$ and $Q$ or tangent to $P$ if $P = Q$.
  \item $\Oinf$ is the neutral element under this law: if 3 points are collinear, their sum is equal to $\Oinf$.
\end{itemize}

This operation are defined in Coq as follow:
\begin{lstlisting}[language=Coq]
Definition neg (p : point) :=
  if p is (| x, y |) then (| x, -y |) else EC_Inf.

Definition add (p1 p2 : point) :=
  match p1, p2 with
    | \infty , _ => p2
    | _ , \infty => p1
    | (| x1, y1 |), (| x2, y2 |) =>
      if x1 == x2 then ... else
        let s := (y2 - y1) / (x2 - x1) in
        let xs := s^2 - x1 - x2 in
          (| xs, - s * (xs - x1 ) - y1 |)
  end.
\end{lstlisting}

And is proven internal to the curve (with coercion):
\begin{lstlisting}[language=Coq]
Lemma addO (p q : ec): oncurve (add p q).

Definition addec (p1 p2 : ec) : ec :=
  EC p1 p2 (addO p1 p2)
\end{lstlisting}

\subsubsection{Montgomery Curves}
\label{montgomery}
Computation over elliptic curves are hard. Speedups can be obtained by using
homogeneous coordinates and other forms than the Weierstra{\ss} form. We consider
the Montgomery form \cite{MontgomerySpeeding}.

\begin{definition}
  Let $a \in \K \backslash \{-2, 2\}$, and $b \in \K \backslash \{ 0\}$.
  The \textit{Montgomery curve} $M_{a,b}(\K)$ is the set of all points $(x,y) \in \K^2$ satisfying the equation:
  $$by^2 = x^3 + ax^2 + x,$$
  along with an additional formal point $\Oinf$, ``at infinity''.
\end{definition}
Using a similar representation, we defined the parametric type \texttt{mc} which
represent the points on a specific Montgomery curve. It is parameterized by
a \texttt{K : ecuFieldType} -- the type of fields which characteristic is not 2 or 3 --
and \texttt{M : mcuType} -- a record that packs the curve parameters $a$ and $b$
along with the proofs that $b \neq 0$ and $a^2 \neq 4$.
\begin{lstlisting}[language=Coq]
Record mcuType :=
  { cA : K; cB : K; _ : cB != 0; _ : cA^2 != 4}.
Definition oncurve (p : point K) :=
if p is (| x, y |)
  then cB * y^+2 == x^+3 + cA * x^+2 + x
  else true.
Inductive mc : Type := MC p of oncurve p.

Lemma oncurve_mc: forall p : mc, oncurve p.
\end{lstlisting}
We define the addition on Montgomery curves the same way as it is in the Weierstra{\ss} form,
however the actual computations will be slightly different.
\begin{lstlisting}[language=Coq]
Definition add (p1 p2 : point K) :=
  match p1, p2 with
    | \infty, _ => p2
    | _, \infty => p1
    | (|x1, y1|), (|x2, y2|) =>
      if   x1 == x2
      then if  (y1 == y2) && (y1 != 0)
           then ... else \infty
      else
        let s  := (y2 - y1) / (x2 - x1) in
        let xs := s^+2 * cB - cA - x1 - x2 in
          (| xs, - s * (xs - x1) - y1 |)
    end.
\end{lstlisting}
But we prove it is internal to the curve (again with coercion):
\begin{lstlisting}[language=Coq]
Lemma addO (p q : mc) : oncurve (add p q).
Definition addmc (p1 p2 : mc) : mc :=
  MC p1 p2 (addO p1 p2)
\end{lstlisting}

We then prove a bijection between a Montgomery curve and its Weierstra{\ss} equation.
\begin{lemma}
  Let $M_{a,b}(\K)$ be a Montgomery curve, define $$a' = \frac{3-a^2}{3b^2} \text{\ \ \ \ and\ \ \ \ } b' = \frac{2a^3 - 9a}{27b^3}.$$
  then $E_{a',b'}(\K)$ is an elliptic curve, and the mapping $\varphi : M_{a,b}(\K) \mapsto E_{a',b'}(\K)$ defined as:
  \begin{align*}
    \varphi(\Oinf_M) &= \Oinf_E\\
    \varphi( (x , y) ) &= ( \frac{x}{b} + \frac{a}{3b} , \frac{y}{b} )
  \end{align*}
  is an isomorphism between groups.
\end{lemma}
\begin{lstlisting}[language=Coq]
Definition ec_of_mc_point p :=
  match p with
  | \infty => \infty
  | (|x, y|) => (|x/(M#b) + (M#a)/(3%:R * (M#b)), y/(M#b)|)
  end.
Lemma ec_of_mc_point_ok p :
  oncurve M p ->
  ec.oncurve E (ec_of_mc_point p).

Definition ec_of_mc p :=
  EC (ec_of_mc_point_ok (oncurve_mc p)).

Lemma ec_of_mc_bij : bijective ec_of_mc.
\end{lstlisting}

\subsubsection{Projective Coordinates}
\label{projective}
Points on a projective plane are represented with a triple $(X:Y:Z)$. Any points
except $(0:0:0)$ defines a point on a projective plane. A scalar multiple of a point defines the same point, \ie
for all $\alpha \neq 0$, $(X:Y:Z)$ and $(\alpha X:\alpha Y:\alpha Z)$ defines
the same point. For $Z\neq 0$, the projective point $(X:Y:Z)$ corresponds to the
point $(X/Z,Y/Z)$ on the Euclidean plane, likewise the point $(X,Y)$ on the Euclidean plane corresponds to $(X:Y:1)$ on the projective plane.

We write the equation for a Montgomery curve $M_{a,b}(\K)$ as such:
\begin{equation}
b \bigg(\frac{Y}{Z}\bigg)^2 = \bigg(\frac{X}{Z}\bigg)^3 + a \bigg(\frac{X}{Z}\bigg)^2 + \bigg(\frac{X}{Z}\bigg)
\end{equation}
Multiplying both sides by $Z^3$ yields:
\begin{equation}
b Y^2Z = X^3 + a X^2Z + XZ^2
\end{equation}
With this equation we can additionally represent the ``point at infinity''. By
setting $Z=0$, we derive $X=0$, giving us the ``infinite points'' $(0:Y:0)$ with $Y\neq 0$.

By restricting the parameter $a$ of $M_{a,b}(\K)$ such that $a^2-4$ is not a
square in \K, we ensure that $(0,0)$ is the only point with a $y$-coordinate of $0$.
\begin{lstlisting}[language=Coq]
Hypothesis mcu_no_square : forall x : K, x^+2 != (M#a)^+2 - 4%:R.
\end{lstlisting}

With those coordinates we prove the following lemmas for the addition of two points.
\begin{definition}Define the functions $\chi$ and $\chi_0$:\\
-- $\chi : M_{a,b}(\K) \to \K \cup \{\infty\}$\\
  such that $\chi(\Oinf) = \infty$ and $\chi((x,y)) = x$.\\
-- $\chi_0 : M_{a,b}(\K) \to \K$\\
  such that $\chi_0(\Oinf) = 0$ and $\chi_0((x,y)) = x$.
\end{definition}
\begin{lemma}
\label{lemma-add}
Let $M_{a,b}(\K)$ be a Montgomery curve such that $a^2-4$ is not a square, and
let $X_1, Z_1, X_2, Z_2, X_3, Z_3 \in \K$, such that $(X_1,Z_1) \neq (0,0)$, $(X_2,Z_2) \neq (0,0)$, $X_4 \neq 0$ and $Z_4 \neq 0$.
Define
\begin{align*}
X_3 &= Z_4((X_1 - Z_1)(X_2+Z_2) + (X_1+Z_1)(X_2-Z_2))^2\\
Z_3 &= X_4((X_1 - Z_1)(X_2+Z_2) - (X_1+Z_1)(X_2-Z_2))^2,
\end{align*}
then for any point $P_1$ and $P_2$ on $M_{a,b}(\K)$ such that $X_1/Z_1 = \chi(P_1), X_2/Z_2 = \chi(P_2)$, and $X_4/Z_4 = \chi(P_1 - P_2)$, we have $X_3/Z_3 = \chi(P_1+P_2)$.\\
\textbf{Remark:} For any $x \in \K \backslash\{0\}, x/0$ should be understood as $\infty$.
\end{lemma}
% This can be formalized as follow:
% \begin{lstlisting}[language=Coq]
% Inductive K_infty :=
% | K_Inf : K_infty
% | K_Fin : K -> K_infty.
%
% Definition point_x (p : point K) :=
%   if p is (|x, _|) then K_Fin x else K_Inf.
% Local Notation "p '#x'" := (point_x p) (at level 30).
% Definition point_x0 (p : point K) :=
%   if p is (|x, _|) then x else 0.
%   Local Notation "p '#x0'" := (point_x0 p) (at level 30).
%
% Definition inf_div (x z : K) :=
%   if z == 0 then K_Inf else K_Fin (x / z).
% Definition hom_ok (x z : K) := (x != 0) || (z != 0).
% Lemma montgomery_hom_neq :
%   forall x1 x2 x4 z1 z2 z4 : K,
%   hom_ok x1 z1 -> hom_ok x2 z2 ->
%   (x4 != 0) && (z4 != 0) ->
%   let x3 := z4 * ((x1 - z1)*(x2 + z2)
%     + (x1 + z1)*(x2 - z2))^+2 in
%   let z3 := x4 * ((x1 - z1)*(x2 + z2)
%     - (x1 + z1)*(x2 - z2))^+2 in
%   forall p1 p2 : point K,
%   oncurve M p1 -> oncurve M p2 ->
%   p1#x = inf_div x1 z1 ->
%   p2#x = inf_div x2 z2 ->
%   (p1 \- p2)#x = inf_div x4 z4 ->
%   hom_ok x3 z3 && ((p1 \+ p2)#x == inf_div x3 z3).
% \end{lstlisting}

With those coordinates we also prove a similar lemma for point doubling.
\begin{lemma}
\label{lemma-double}
Let $M_{a,b}(\K)$ be a Montgomery curve such that $a^2-4$ is not a square, and
let $X_1, Z_1, X_2, Z_2, X_3, Z_3 \in \K$, such that $(X_1,Z_1) \neq (0,0)$. Define
\begin{align*}
  c &= (X_1 + Z_1)^2 - (X_1 - Z_1)^2\\
X_3 &= (X_1 + Z_1)^2(X_1-Z_1)^2\\
Z_3 &= c\Big((X_1 + Z_1)^2+\frac{a-2}{4}\times c\Big),
\end{align*}
then for any point $P_1$ on $M_{a,b}(\K)$ such that $X_1/Z_1 = \chi(P_1)$,
we have $X_3/Z_3 = \chi(2P_1)$.
\end{lemma}
% Which is formalized as follow:
% \begin{lstlisting}[language=Coq]
% Lemma montgomery_hom_eq :
%   forall x1 z1 : K,
%   hom_ok x1 z1 ->
%   let c := (x1 + z1)^+2 - (x1 - z1)^+2 in
%   let x3 := (x1 + z1)^+2 * (x1 - z1)^+2 in
%   let z3 := c * ((x1 + z1)^+2 + (((M#a) - 2%:R)/4%:R) * c) in
%   forall p : point K, oncurve M p ->
%   p#x = inf_div x1 z1 ->
%   (p \+ p)#x = inf_div x3 z3.
% \end{lstlisting}

With these two lemmas (\ref{lemma-add} and \ref{lemma-double}), we have the basic
tools to compute efficiently additions and point doubling on projective coordinates.

\subsubsection{Scalar Multiplication Algorithms}
\label{ladder}

Suppose we have a scalar $n$ and a point $P$ on some curve. The most straightforward
way to compute $n \cdot P$ is to repetitively add $P$ \ie computing $P + \ldots + P$.
However there is an more efficient algorithm which makes use of the binary
representation of $n$ and by combining doubling and adding and starting from $\Oinf$.
\eg for $n=11$, we compute $2(2(2(2\Oinf + P)) + P)+ P$.

% \begin{algorithm}
% \caption{Double-and-add for scalar mult.}
% \label{double-add}
% \begin{algorithmic}
% \REQUIRE{Point $P$, scalars $n$ and $m$, $n < 2^m$}
% \ENSURE{$Q = n \cdot P$}
% \STATE $Q \leftarrow \Oinf$
% \FOR{$k$ := $m$ downto $1$}
%   \STATE $Q \leftarrow 2Q$
%   \IF{$k^{\text{th}}$ bit of $n$ is $1$}
%     \STATE $Q \leftarrow Q + P$
%   \ENDIF
% \ENDFOR
% \end{algorithmic}
% \end{algorithm}

% \begin{lemma}
% \label{lemma-double-add}
% Algorithm \ref{double-add} is correct, \ie it respects its output conditions given the input conditions.
% \end{lemma}

% We prove Lemma \ref{lemma-double-add}. However
With a simple double-and-add algorithm, with careful timing, an attacker could reconstruct $n$.
In the case of X25519, $n$ is the private key. With the Montgomery's ladder, while
it provides slightly more computations and an extra variable, we can prevent such weakness.
See Algorithm \ref{montgomery-ladder}.

% \begin{algorithm}
% \caption{Montgomery ladder for scalar mult.}
% \label{montgomery-ladder}
% \begin{algorithmic}
% \REQUIRE{Point $P$, scalars $n$ and $m$, $n < 2^m$}
% \ENSURE{$Q = n \cdot P$}
% \STATE $Q \leftarrow \Oinf$
% \STATE $R \leftarrow P$
% \FOR{$k$ := $m$ downto $1$}
%   \IF{$k^{\text{th}}$ bit of $n$ is $0$}
%     \STATE $R \leftarrow Q + R$
%     \STATE $Q \leftarrow 2Q$
%   \ELSE
%     \STATE $Q \leftarrow Q + R$
%     \STATE $R \leftarrow 2R$
%   \ENDIF
% \ENDFOR
% \end{algorithmic}
% \end{algorithm}

\begin{lemma}
\label{lemma-montgomery-ladder}
Algorithm \ref{montgomery-ladder} is correct, \ie it respects its output conditions given the input conditions.
\end{lemma}

In Curve25519 we are only interested in the $x$ coordinate of points, using
Lemmas \ref{lemma-add} and \ref{lemma-double}, and replacing the if statements
with conditional swapping we can define a ladder similar to the one used in TweetNaCl.
See Algorithm \ref{montgomery-double-add}

\begin{algorithm}
\caption{Montgomery ladder for scalar multiplication on $M_{a,b}(\K)$ with optimizations}
\label{montgomery-double-add}
\begin{algorithmic}
\REQUIRE{$x \in \K\backslash \{0\}$, scalars $n$ and $m$, $n < 2^m$}
\ENSURE{$a/c = \chi_0(n \cdot P)$ for any $P$ such that $\chi_0(P) = x$}
\STATE $(a,b,c,d) \leftarrow (1,x,0,1)$
\FOR{$k$ := $m$ downto $1$}
  \IF{$k^{\text{th}}$ bit of $n$ is $1$}
    \STATE $(a,b) \leftarrow (b,a)$
    \STATE $(c,d) \leftarrow (d,c)$
  \ENDIF
  \STATE $e \leftarrow a + c$
  \STATE $a \leftarrow a - c$
  \STATE $c \leftarrow b + d$
  \STATE $b \leftarrow b - d$
  \STATE $d \leftarrow e^2$
  \STATE $f \leftarrow a^2$
  \STATE $a \leftarrow c \times a$
  \STATE $c \leftarrow b \times e$
  \STATE $e \leftarrow a + c$
  \STATE $a \leftarrow a - c$
  \STATE $b \leftarrow a^2$
  \STATE $c \leftarrow d-f$
  \STATE $a \leftarrow c\times\frac{A - 2}{4}$
  \STATE $a \leftarrow a + d$
  \STATE $c \leftarrow c \times a$
  \STATE $a \leftarrow d \times f$
  \STATE $d \leftarrow b \times x$
  \STATE $b \leftarrow e^2$
  \IF{$k^{\text{th}}$ bit of $n$ is $1$}
    \STATE $(a,b) \leftarrow (b,a)$
    \STATE $(c,d) \leftarrow (d,c)$
  \ENDIF
\ENDFOR
\end{algorithmic}
\end{algorithm}

\begin{lemma}
\label{lemma-montgomery-double-add}
Algorithm \ref{montgomery-double-add} is correct, \ie it respects its output
conditions given the input conditions.
\end{lemma}

%% here we have \chi and \chi_0 ...

We formalized this lemma (\ref{lemma-montgomery-double-add}):
\begin{lstlisting}[language=Coq]
Lemma opt_montgomery_x :
  forall (n m : nat) (x : K),
  n < 2^m -> x != 0 ->
  forall (p : mc M), p#x0 = x ->
  opt_montgomery n m x = (p *+ n)#x0.
\end{lstlisting}

We can remark that for an input $x = 0$, the ladder returns $0$.
\begin{lstlisting}[language=Coq]
Lemma opt_montgomery_0:
  forall (n m : nat), opt_montgomery n m 0 = 0.
\end{lstlisting}
Also \Oinf\ is the neutral element over $M_{a,b}(\K)$, we have:
$$\forall P, P + \Oinf\ = P$$
thus we derive the following lemma.
% \begin{lemma}
% \label{lemma-montgomery-double-add}
% Algorithm \ref{montgomery-double-add} is correct even if $x=0$, \ie it respects its output conditions given the input conditions or $x=0$.
% \end{lemma}
\begin{lstlisting}[language=Coq]
Lemma p_x0_0_eq_0 : forall (n : nat) (p : mc M),
  p #x0 = 0%:R -> (p *+ n) #x0 = 0%R.
\end{lstlisting}
And thus the theorem of the correctness of the Montgomery ladder.
\begin{theorem}
\label{montgomery-ladder-correct}
For all $n, m \in \N$, $x \in \K$, $P \in M_{a,b}(\K)$,
if $\chi_0(P) = x$ then Algorithm \ref{montgomery-double-add} returns $\chi_0(n \cdot P)$
\end{theorem}
\begin{lstlisting}[language=Coq]
Theorem opt_montgomery_ok (n m: nat) (x : K) :
  n < 2^m ->
  forall (p : mc M), p#x0 = x ->
  opt_montgomery n m x = (p *+ n)#x0.
\end{lstlisting}

\subsection{Curves, Twists and Extension Fields}

One hypothesis to be able to use the above theorem is that $a^2-4$ is not a square:
$$\forall x \in \K,\ x^2 \neq a^2-4$$
As Curve25519 is defined over the field $\K = \F{p^2}$, there exists $x$ such that $x^2 = a^2-4$.
We first study Curve25519 and one of the quadratic twist Twist25519, first defined over \F{p}.

\subsubsection{Curves and Twists}

We define $\F{p}$ as the numbers between $0$ and $p = \p$.
We create a \coqe{Zmodp} module to encapsulate those definitions.
\begin{lstlisting}[language=Coq]
Module Zmodp.
Definition betweenb x y z := (x <=? z) && (z <? y).
Definition p := locked (2^255 - 19).
Fact Hp_gt0 : p > 0.
Inductive type := Zmodp x of betweenb 0 p x.

Lemma Z_mod_betweenb (x y : Z) :
  y > 0 -> betweenb 0 y (x mod y).

Definition pi (x : Z) : type :=
  Zmodp (Z_mod_betweenb x Hp_gt0).
Coercion repr (x : type) : Z :=
  let: @Zmodp x _ := x in x.
End Zmodp.
\end{lstlisting}

We define the basic operations ($+, -, \times$) with their respective neutral elements ($0, 1$).
\begin{lemma}
$\F{p}$ is a commutative ring.
\end{lemma}
% \begin{lstlisting}[language=Coq]
% Definition zero : type := pi 0.
% Definition one : type := pi 1.
% Definition opp (x : type) : type := pi (p - x).
% Definition add (x y : type) : type := pi (x + y).
% Definition sub (x y : type) : type := pi (x - y).
% Definition mul (x y : type) : type := pi (x * y).
%
% Lemma Zmodp_ring :
%   ring_theory zero one add mul sub opp eq.
% \end{lstlisting}
And finally for $a = 486662$, by using the Legendre symbol we prove that $a^2 - 4$ and $2$ are not squares in $\F{p}$.
\begin{lstlisting}[language=Coq]
Lemma a_not_square : forall x: Zmodp.type,
  x^+2 != (Zmodp.pi 486662)^+2 - 4%:R.
\end{lstlisting}
\begin{lstlisting}[language=Coq,label=two_not_square]
Lemma two_not_square : forall x : Zmodp.type,
  x^+2 != 2%:R.
\end{lstlisting}
We consider $M_{486662,1}(\F{p})$ and $M_{486662,2}(\F{p})$, one of its quadratic twist.
% $M_{486662,1}(\F{p})$ has the same equation as $M_{486662,1}(\F{p^2})$ while $M_{486662,2}(\F{p})$ is one of its quadratic twist.


By instantiating theorem \ref{montgomery-ladder-correct} we derive the following two lemmas:
\begin{lemma} For all $x \in \F{p},\ n \in \N,\ P \in \F{p} \times \F{p}$,\\
such that $P \in M_{486662,1}(\F{p})$ and $\chi_0(P) = x$.
Given $n$ and $x$, $Curve25519\_Fp(n,x) = \chi_0(n \cdot P)$.
\end{lemma}
\begin{lemma} For all $x \in \F{p},\ n \in \N,\ P \in \F{p} \times \F{p}$\\
such that $P \in M_{486662,2}(\F{p})$ and $\chi_0(P) = x$.
Given $n$ and $x$, $Twist25519\_Fp(n,x) = \chi_0(n \cdot P)$.
\end{lemma}
As the Montgomery ladder defined above does not depends on $b$, it is trivial to see that the computations done for points of $M_{486662,1}(\F{p})$ and of $M_{486662,2}(\F{p})$ are the same.
\begin{lstlisting}[language=Coq]
Theorem curve_twist_eq: forall n x,
  curve25519_Fp_ladder n x = twist25519_Fp_ladder n x.
\end{lstlisting}

Because $2$ is not a square in $\F{p}$, it allows us split $\F{p}$ into two sets.
\begin{lemma}
\label{square-or-2square}
For all $x$ in $\F{p}$, there exists $y$ in $\F{p}$ such that
$$y^2 = x\ \ \ \lor\ \ 2y^2 = x$$
\end{lemma}
For all $x \in \F{p}$, we can compute $x^3 + ax^2 + x$. Using Lemma \ref{square-or-2square} we can find a $y$ such that $(x,y)$ is either on the curve or on the quadratic twist:
\begin{lemma}
\label{curve-or-twist}
For all $x \in \F{p}$, there exists a point $P$ over $M_{486662,1}(\F{p})$ or over $M_{486662,2}(\F{p})$ such that the $x$-coordinate of $P$ is $x$.
\end{lemma}
\begin{lstlisting}[language=Coq]
Theorem x_is_on_curve_or_twist: forall x : Zmodp.type,
  (exists (p : mc curve25519_mcuType), p#x0 = x) \/
  (exists (p' : mc twist25519_mcuType), p'#x0 = x).
\end{lstlisting}

\subsubsection{Curve25519 over \F{p^2}}

We use the same definitions as in \cite{Ber06}. We consider the extension field $\F{p^2}$ as the set $\F{p} \times \F{p}$ with $\delta = 2$, in other words,
the polynomial with coefficients in $\F{p}$ modulo $X^2 - 2$. In a similar way as for $\F{p}$ we use Module in Coq.
\begin{lstlisting}[language=Coq]
Module Zmodp.
Inductive type :=
  Zmodp2 (x: Zmodp.type) (y:Zmodp.type).

Definition pi (x : Zmodp.type * Zmodp.type) : type :=
  Zmodp2 x.1 x.2.
Coercion repr (x : type) : Zmodp.type*Zmodp.type :=
  let: Zmodp2 u v := x in (u, v).

Definition zero : type :=
  pi ( 0%:R, 0%:R ).
Definition one : type :=
  pi ( 1, 0%:R ).
Definition opp (x : type) : type :=
  pi (- x.1 , - x.2).
Definition add (x y : type) : type :=
  pi (x.1 + y.1, x.2 + y.2).
Definition sub (x y : type) : type :=
  pi (x.1 - y.1, x.2 - y.2).
Definition mul (x y : type) : type :=
  pi ((x.1 * y.1) + (2%:R * (x.2 * y.2)),
      (x.1 * y.2) + (x.2 * y.1)).
\end{lstlisting}
We define the basic operations ($+, -, \times$) with their respective neutral elements ($0, 1$).
Additionally we verify that for each element of in $\F{p^2}\backslash\{0\}$, there exists a multiplicative inverse.
\begin{lemma} For all $x \in \F{p^2}\backslash\{0\}$ and $a,b \in \F{p}$ such that $x = (a,b)$,
$$x^{-1} = \Big(\frac{a}{a^2-2b^2}\ , \frac{-b}{a^2-2b^2}\Big)$$
\end{lemma}
Similarly as in $\F{p}$, we define $0^{-1} = 0$.
\begin{lemma}
$\F{p^2}$ is a commutative ring.
\end{lemma}
We can then specialize the basic operations in order to speed up the verification of formulas by using rewrite rules:
\begin{align*}
(a,0) + (b,0) &= (a+b, 0)\\
(a,0) \cdot   (b,0) &= (a \cdot b, 0)\\
(a, 0)^n &= (a^n, 0)\\
(a, 0)^{-1} &= (a^{-1}, 0)\\
(a, 0)\cdot (0,b) &= (0, a\cdot b)\\
(0, a)\cdot (0,b) &= (2\cdot a\cdot b, 0)\\
(0,a)^{-1} &= ((2\cdot a)^{-1},0)
\end{align*}
The injection $a \mapsto (a,0)$ from $\F{p}$ to $\F{p^2}$ preserves $0, 1, +, -, \times$. Thus $(a,0)$ can be abbreviated as $a$ without confusions.

We define $M_{486662,1}(\F{p^2})$. With the rewrite rule above, it is straightforward to prove that any point on the curve $M_{486662,1}(\F{p})$ is also on the curve $M_{486662,1}(\F{p^2})$. Similarly, any point on the quadratic twist $M_{486662,2}(\F{p})$ is also on the curve $M_{486662,1}(\F{p^2})$.
As direct consequence, using lemma \ref{curve-or-twist}, we prove that for all $x \in \F{p}$, there exists a point $P \in \F{p^2}\times\F{p^2}$ on $M_{486662,2}(\F{p})$ such that $\chi_0(P)$ is $(x,0)$

\begin{lstlisting}[language=Coq]
Theorem x_is_on_curve_or_twist_implies_x_in_Fp2:
  forall (x:Zmodp.type),
    exists (p: mc curve25519_Fp2_mcuType),
      p#x0 = Zmodp2.Zmodp2 x 0.
\end{lstlisting}

We now study the case of the scalar multiplication and show similar proofs.
\begin{definition}
Define the functions $\varphi_c$, $\varphi_t$ and $\psi$\\
-- $\varphi_c: M_{486662,1}(\F{p}) \mapsto M_{486662,1}(\F{p^2})$\\
  such that $\varphi((x,y)) = ((x,0), (y,0))$.\\
-- $\varphi_t: M_{486662,2}(\F{p}) \mapsto M_{486662,1}(\F{p^2})$\\
  such that $\varphi((x,y)) = ((x,0), (0,y))$.\\
-- $\psi: \F{p^2} \mapsto \F{p}$\\
  such that $\psi(x,y) = (x)$.
\end{definition}

\begin{lemma}
For all $n \in \N$, for all point $P\in\F{p}\times\F{p}$ on the curve $M_{486662,1}(\F{p})$ (respectively on the quadratic twist $M_{486662,2}(\F{p})$), we have:
\begin{align*}
P \in M_{486662,1}(\F{p}) &\implies \varphi_c(n \cdot P) = n \cdot \varphi_c(P)\\
P \in M_{486662,2}(\F{p}) &\implies \varphi_t(n \cdot P) = n \cdot \varphi_t(P)
\end{align*}
\end{lemma}

Notice that:
\begin{align*}
\forall P \in M_{486662,1}(\F{p}),\ \ \psi(\chi_0(\varphi_c(P))) = \chi_0(P)\\
\forall P \in M_{486662,2}(\F{p}),\ \ \psi(\chi_0(\varphi_t(P))) = \chi_0(P)
\end{align*}

In summary for all $n \in \N,\ n < 2^{255}$, for any given point $P\in\F{p}\times\F{p}$ on $M_{486662,1}(\F{p})$ or $M_{486662,2}(\F{p})$ \texttt{curve25519\_Fp\_ladder} computes the $\chi_0(n \cdot P)$.
We have proved that for all $P \in \F{p^2}\times\F{p^2}$ such that $\chi_0(P) \in \F{p}$ there exists a corresponding point on the curve or the twist over $\F{p}$.
We have proved that for any point, on the curve or the twist we can compute the scalar multiplication by $n$ and yield to the same result as if we did the computation in $\F{p^2}$. As a result we have proved theorem 2.1 of \cite{Ber06}:
\begin{theorem}
For all $n \in \N$, $x \in \F{P}$, $P \in M_{486662,1}(\F{p^2})$, such that $n < 2^{255}$ and $\chi_0(P) = \varphi(x)$, \texttt{curve25519\_Fp\_ladder}$(n, x)$ computes $\psi(\chi_0(n \cdot P))$.
\end{theorem}
which can be formalized in Coq as:
\begin{lstlisting}[language=Coq]
Lemma curve25519_Fp2_ladder_ok (n : nat) (x:Zmodp.type) :
    (n < 2^255)%nat ->
    forall (p  : mc curve25519_Fp2_mcuType),
    p #x0 = Zmodp2.Zmodp2 x 0 ->
    curve25519_Fp_ladder n x = (p *+ n)#x0 /p.
\end{lstlisting}