**Input:** a biextension element $g_{P,Q} = [\tilde{P}, \tilde{Q}, \tilde{\theta}, P + Q]$ represented by cubical points

**Output:** cubical points $\ell\tilde{P}, \ell\tilde{P} + Q$ such that $g_{P,Q}^{*1,\ell} = [\ell\tilde{P}, \tilde{Q}, \tilde{\theta}, \ell\tilde{P} + Q]$

→ For each bit $b_i$ of $\ell$ from left to right (skipping the first one), given a cubical representation $\widetilde{nP}, n\widetilde{P} + Q$, of $g_{P,Q}^{*1,n}$, where $n$ is the current truncation of $\ell$ on the leftmost bits, do a cubical biextension double:

a. $2n\widetilde{P} + Q = \text{DiffAdd}(n\widetilde{P} + Q, n\widetilde{P}, \tilde{Q})$;

b. $2n\widetilde{P} = \text{Double}(n\widetilde{P})$;

And if $b_i = 1$, also do a cubical biextension addition:

a. Compute $(2n+1)P = \text{CompatibleAdd}(2nP, P; 2nP+Q, P+Q)$ and take an arbitrary cubical lift $(2n+1)P$ of $(2n+1)P$.

b. $(2n+1)P + Q = \text{ThreeWayAdd}(2n\widetilde{P}, \tilde{P}, \tilde{Q}, P + Q, 2n\widetilde{P} + Q, (2n+1)P)$;

<caption>ALGORITHM 4.7. Biextension exponentiation on Kummer varieties via cubical double and add</caption>

$X_i(P \pm Q)$ and $X_1(P \pm Q)$. One can then consider the abelian and Kummer schemes $A_R/R$, $(A_R / \pm 1)/R$ over $R$. The coordinate $X \bmod \mathfrak{P}$ encodes $P + Q$ and $P - Q$. More formally, we have a canonical $R$ point Spec $R \to A_R / \pm 1$, and as, eventually over a degree 2 extension $k'$ of $k$, $R$ splits as $R_{k'} = k' \oplus k'$, the isomorphism Spec $R \cong \pm(P \pm Q)$ splits as Spec $k' \oplus$ Spec $k' \cong \pm(P + Q) \times \pm(P - Q)$, so the canonical $R$ point splits as the point $\pm(P + Q)$ and the point $\pm(P - Q)$. Anyway, working over $R$ rather than over $k$, we can do our pairing computations as before, except that in the end we obtain a representation of the degree two scheme $e(P, Q)^{\pm 1}$ of $\mathbb{G}_m$. Typically this representation is given by a degree two polynomial $\mathfrak{Q}(X) = (X - e(P, Q))(X - e(P, Q)^{-1}) = X^2 - (e(P, Q) + e(P, Q)^{-1})X + 1$, and so we can recover the trace $(e(P, Q) + e(P, Q)^{-1})$, which we will call the symmetric pairing of $P, Q$.

Alternatively we can see $e(P, Q)^{\pm 1}$ as a point of $\mathbb{G}_m / \pm 1$, and if $x^{\pm 1} \in \mathbb{G}_m / \pm 1$, the trace $x^{\pm 1} \mapsto x + 1/x$ is a convenient representation of $x^{\pm 1}$. We can still do arithmetic on $\mathbb{G}_m / \pm 1$ (see [Rob21a, § 2.12.2]), notably compute exponentiation via squarings and differential multiplications. This allows us to compute reduced symmetric pairings from non reduced symmetric pairings. This strategy to compute pairings is well known, see [GL08] for Kummer lines and [LR15] for the case of Kummer varieties in the theta model.

### 4.8. Analytic cubical points and analytic theta functions.
Our goal is now to find formulas for the cubical arithmetic using the affine lift representation. It will be convenient to work out formulas over $\mathbb{C}$.

Since we know that the cubical arithmetic is algebraic, working over the universal abelian scheme there certainly exists algebraic formulas over $\mathbb{C}$. By standard arguments (see [Rob21b, § 2.3.6]), if we find through analytic means algebraic formulas defiend over $\mathbb{Q}$, we know that they give the correct arithmetic formula over any field $k$ of characteristic $p$, as long as the formulas have good reduction modulo $p$.

Now, from the algebraic Riemann relations on line bundles from Proposition 4.1, taking sections $s_i$ of $\mathcal{L}$ we know that there should exists linear realations between suitable translated tensor products of the $s_i$, i.e. Riemann relations on sections. These have been worked out in the case where the $s_i$ are the algebraic theta functions by Mumford in [Mum66, p. 333–335].

Now let $A/\mathbb{C}$ be an abelian variety, it is known that it is a (polarisable) torus $A = V/\Lambda$ with $V$ a $\mathbb{C}$-vector space of dimension $g$, and that the analytic addition law on $V \simeq \mathbb{C}^g$
induces the algebraic addition law on $\mathcal{A}$. Fix a very ample line bundle $\mathcal{L}$ on $\mathcal{A}$. It is also
known that there is a theory of analytic theta functions on $\mathcal{A}$ with respect to $\mathcal{L}$, which are
given by analytic functions on $V$ [Mum83; Mum84].

Given these theta functions $\theta_i$, we have two representations of a point $P \in \mathcal{A}$. First we have the projective representation: take any $z_P \in V$ above $P$, and let $\theta(P) = (\theta_1(z_P) : \dots : \theta_m(z_P))$. The theta functions are automorphic with respect to $\Lambda$ (with the same factor of automorphy), so changing $z_P$ does not change the projective point $\theta(P)$. But we also have an affine representation $\theta(z_P) = (\theta_1(z_P), \dots, \theta_m(z_P))$, which is an affine point above $\theta(P)$.

We remark the similarity with the discussion of Section 4.5. In fact, since $V$ is simply connected, the pullback of any line bundle $\mathcal{L}$ on $\mathcal{A}$ to $V$ is trivial, so essentially the choice of $z_P$ above $P$ induces a choice of rigidification of $\mathcal{L}$ at $P$.

If we know the projective points $\theta(P), \theta(Q)$, we can compute $\theta(P + Q)$. However, if we know $\theta(z_P), \theta(z_Q)$, we cannot recover the analytic addition $\theta(z_P + z_Q)$. However, we can use the analytic Riemann relations [Mum83], so that whenever we have $[z_1, z_2, z_3, z_4; z'_1, z'_2, z'_3, z'_4]$ in Riemann position, and we know all but one of the $\theta(z_i)$ we can recover exactly the last one. In particular, as in Example 4.3, we can do analytic differential additions and analytic three way additions.

Now the analytic Riemann relations on analytic theta functions are exactly the same as the algebraic Riemann relations on algebraic theta functions. Since the analytic Riemann relations are induced by the analytic group law on $V$ on one hand, and the algebraic Riemann relations give the cubical structure, we deduce that the analytic addition law on $V$ induces the cubical algebraic structure. See also [Bre83, p. 41, 42]. In fact Breen uses the algebraic cubical theory to give an alternative definition of algebraic theta functions compared to Mumford (which uses the Heisenberg group representation and a choice of isomorphism between the Heisenberg group and a theta group).

In summary, the cubical arithmetic encodes the algebraic information which can be extracted from the analytic group law of $V$ above the algebraic group law of $\mathcal{A}$. Analytically, a cubical point $\bar{P}$ corresponds to an analytic point $z_P$ above $P$, and cubical arithmetic like $m\bar{P} + nQ$ (when we have enough information to compute it) corresponds to the analytic point $mz_P + nz_Q$. This allows us to naturally check the various compatibility relations on the cubical arithmetic by checking it through the analytic group law.

The first consequence, is that the affine version of the theta Riemann relations, since they are over $\mathbb{C}$ given by the analytic theta Riemann relations coming from the analytic group law, give explicit formulas for the cubical arithmetic expressed in terms of affine lifts of theta points. In particular, all the arithmetic on affine lift of theta points as developed in [LR16; LR10; LR15], was actually a cubical arithmetic in disguise, as should be clear when comparing Theorem 4.19 with [LR16; LR15]. This explains how to do the cubical arithmetic in theta models of abelian varieties. This allows to have a nice intuition on the operations allowed under the cubical arithmetic: any algebraic relations that can be derived from the transcendental addition law $z_i + z_j + \cdots + z_k$ via the analytic theta functions (so for instance: we cannot recover the $\theta_i(z_1 + z_2)$ from the $\theta_i(z_1)$, $\theta_i(z_2)$, but we can is we also know the $\theta_i(z_1 - z_2)$ can be rewritten as some cubical operation.

Furthermore, the automorphic factor associated to theta functions gives an explicit trivialisation of the theta groups on $A = V/\Lambda$ pulled back to $V$, and of biextensions on $A \times A$ pulled back to $V \times V$. This allows to give an analytic proof of the link between cubical points and biextensions.

We conclude this section by an example of how the analytic interpretation of the cubical
law makes it very easy to prove certain of its properties:
**Lemma 4.32.** Let $M$ be an $8 \times 8$ matrix of cubical points. Assume that all lines give cubical points in Riemann position, and 7 out of 8 columns also give cubical points in Riemann position. Then the last column also gives cubical points in Riemann position.

*Proof.* Let $[\tilde{P}_1, \tilde{P}_2, \tilde{P}_3, \tilde{P}_4; \tilde{Q}_1, \tilde{Q}_2, \tilde{Q}_3, \tilde{Q}_4]$ be analytic cubical points in Riemann position. This means that we have an affine representation via analytic coordinates $\tilde{P}_i = (\theta_j(z_i))$, $\tilde{Q}_i = (\theta_j(z'_i))$, where $z_i, z'_i \in \mathbb{C}^g$ and $[z_1, z_2, z_3, z_4; z'_1, z'_2, z'_3, z'_4]$ are analytic points in Riemann position.

To our matrix $M$ of cubical points, corresponds a matrix $M_z$ of elements $z_{ij} \in \mathbb{C}^g$, such that each line and 7 out of 8 columns are in Riemann position. Then by linear algebra the last column has to be in Riemann position too.

In other words: going from the $\tilde{P}_i$ to the analytic $z_i$ is a way to trivialize the cubical arithmetic; for an algebraic proof, we can use [Mor85, § I.5] instead. $\square$

**Corollary 4.33.** If $[\tilde{P}_1, \tilde{P}_2, \tilde{P}_3, \tilde{P}_4; \tilde{Q}_1, \tilde{Q}_2, \tilde{Q}_3, \tilde{Q}_4]$ are cubical points in Riemann position, so are $[m\tilde{P}_1, m\tilde{P}_2, m\tilde{P}_3, m\tilde{P}_4; m\tilde{Q}_1, m\tilde{Q}_2, m\tilde{Q}_3, m\tilde{Q}_4]$ for any $m \in \mathbb{Z}$.

**4.9. Cubical arithmetic on elliptic curves and Kummer lines.** In this section, we give explicit formulas for the arithmetic of cubical elliptic points in Sections 4.9.1 to 4.9.3, and for cubical points on elliptic Kummer lines in Section 4.9.4. Finally in Section 4.9.5 we make the link between elliptic nets and cubical points.

**4.9.1. Weierstrass coordinates.** We will use Section 4.8 to derive the cubical arithmetic of an elliptic curve $(E, (0_E))$ in Weierstrass equation $y^2 = x^3 + a_2x^2 + a_4x + a_6$.

It is well known by Riemann Roch that the line bundle $\mathcal{L}$ associated to $(0_E)$ has one global section $Z_1$. The line bundle $\mathcal{L}^2$ has two sections, $X_2, Z_2$ with $Z_2 = Z_1^2$. And finally the line bundle $\mathcal{L}^3$ has three sections $X_3, Y_3, Z_3$, with $X_3 = X_2Z_1, Z = Z_1^3$.

To simplify notations, we let $X, Y, Z = X_3, Y_3, Z_3$. We can work on the elliptic curve via the projective coordinates $(X : Y : Z)$, or via the affine coordinates $x = X/Z = X_2/Z_2, y = Y/Z$ which are defined everywhere except at $0_E$. (It is also standard in elliptic curve cryptography to use Jacobian coordinates which are given by $(X_2 : Y_3 : Z_1)$, with projective weights $(2 : 3 : 1)$.)

Since $0_E$ is of multiplicity 1 in the divisor $(0_E)$, by the proof of [Mum66, Proposition 2 p.307], we have that $Z_1(-\tilde{P}) = -Z_1(\tilde{P})$. So $Z_1$ is odd, and since $x$ is odd and $y$ even, $X$ is even and $Y$ is even.

The affine representation of a cubical point $\tilde{P}$ above $P = (X(P) : Y(P) : Z(P))$, with respect to the divisor $D = 3(0_E)$, then corresponds to a choice of affine lift $\tilde{P} = (X(\tilde{P}), Y(\tilde{P}), Z(\tilde{P}))$ above the projective point $P$.

If we want to work with the cubical arithmetic associated to $(0_E)$ (hence a principal polarisation), since it is not very ample, to represent a cubical point $\tilde{P}$, we need to fix both a representation of $P$ itself, e.g. via the Weierstrass coordinates $(X(P) : Y(P) : Z(P))$ and a choice of $Z_1(P)$, which we will denote by $Z_1(\tilde{P})$. Our representation is thus $\tilde{P} = (P, Z_1(\tilde{P}))$. It is well defined except when $P = 0_E$ because the neutral point is a base point of $(0_E)$ and $Z_1(0_E) = 0$.

**Remark 4.34.** As discussed in Remark 4.24, we want to normalize $\tilde{0}$ such that $(Z_1/(x/y))(\tilde{0}) = 1$. With this normalisation, we have $X_2(\tilde{0}) = (X_2/(Z_1/(x/y))^2)(\tilde{0}) = (X_2x^2/(Z_1^2y^2))(\tilde{0}) = (x^3/y^2)(\tilde{0}) = 1$, and $Y_3(\tilde{0}) = (Y_3/(Z_1/(x/y))^3)(\tilde{0}) = (Y_3x^3/(Z_1^3y^3))(\tilde{0}) = (x^3/y^2)(\tilde{0}) = 1$.
We will call sections of the divisor $n(0_E)$ coordinates of level $n$, so for instance $Z_1$ is of level 1, $X_2$, $Z_2$ of level 2, and $X_3$, $Y_3$, $Z_3$ of level 3. To specify a cubical point $\tilde{P}$ with respect to the divisor $n(0_E)$, we need to specify $P$ using some projective coordinates of level $m$: $(X_{m,i}(P))$, and one or several affine coordinates of level $n$: $X_{n,i}(\tilde{P})$. If $m=n$, then the projective coordinates are subsumed by the affine coordinates, but we can take $m$ different from $n$.

For instance the representation $\tilde{P} = (X(\tilde{P}), Y(\tilde{P}), Z(\tilde{P}))$ uses level 3 affine coordinates (so gives cubical points associated to $3(0_E)$); and our normalisation of the neutral point is $\tilde{0} = (0,1,0)$. The representation $\tilde{P} = ((X(P) : Y(P) : Z(P)), Z_1(\tilde{P}))$ uses level 3 projective coordinates with a level 1 affine coordinate (so cubical points associated to $(0_E)$). This times, $Z_1(0_E) = 0$ so to define our neutral point we need to add the extra condition that $(Z_1/(x/y))(\tilde{0}) = 1$.

In Section 4.9.4 we will use level 2 affine coordinates $\tilde{P} = (X_2(\tilde{P}), Z_2(\tilde{P}))$. Projectively, these only allow to recover $(X_2(P) : Z_2(P))$, i.e. $x(P)$, so they cannot distinguish between $P$ and $-P$. Our neutral point is $\tilde{0} = (1,0)$.

On a twisted Edwards curve, the completed Edward coordinates are given by $\{(X : Z), (Y : T) \in \mathbb{P}^1 \times \mathbb{P}^1 | aX^2T^2 + Y^2Z^2 = Z^2T^2 + dX^2Y^2\}$, and the Segre embedding $\mathbb{P}^1 \times \mathbb{P}^1 \to \mathbb{P}^3$, $(X : Z), (Y : T) \mapsto (XT : YZ : ZT : XY)$ gives the extended Edwards coordinates (which are of level 4) [BBLP13, § 2.7]. Let $M$ be the Montgomery model which is birationally equivalent to $E$, and $X_M, Z_M$ its level two coordinates, and $T = (1:0)$ the canonical point of 4 torsion on the Kummer line. Then up to a linear change of variable, the completed Edwards coordinates correspond precisely to the embedding $M \to \mathbb{P}^1 \times \mathbb{P}^1, P \mapsto ((X_M(P) : Z_M(P)), (X_M(P+T) : Z_M(P+T)))$. This was first remarked (implicitly) in [Koh11b; Koh11a, §8.1]; see also [FK22, § 3] and [LR16] for generalisations to higher dimension. This means we could use a mix of affine level 2 and projective level 4 coordinates to describe cubical points of level 2 in $\mathbb{A}^2 \times \mathbb{P}^1$ via $\tilde{P} = ((X_M(\tilde{P}), Z_M(\tilde{P})), (X_M(P+T) : Z_M(P+T)))$, and compute pairings using a mix of Kummer line cubical arithmetic as in Section 5 and Edwards additions. We leave that for future work.

More generally, on an abelian variety $A$, with projective embedding $A \to \mathbb{P}^N$ given by a very ample divisor $E$, and a divisor $D$, we just need an explicit version of the theorem of the square for $D$, i.e., for $P, Q \in A$, represented by their projective coordinates $X_i(P), X_i(Q)$, $X_i \in \Gamma(E)$, build a function $g_{P,Q}$ with divisor $D_{P+Q} - D_P - D_Q$, such that we are able to evaluate it on $R = (X_i(R))$, to be able to work out the cubical arithmetic, hence biextension arithmetic, associated to $D$, while working in the model associated to $E$. It is not required to choose $E = D$.

<footnote n="4">4.9.2. Analytic cubical elliptic points. Analytically, the Weierstrass sigma function $\sigma$ is a theta function (up to some exponential factor) which has zeroes of order 1 exactly at the lattice $\Lambda$, it can thus play the role of our $Z_1$ (see Remark 4.35). It is also well known that analytically, $x = \wp, y = \wp'$.</footnote>
The cubical arithmetic with respect to $\sigma$ is then given by the Frobenius Stickelberger relations (see [Stao8, Lemma 5.1.3; Sil86, Exercise 6.3; Bre83, Eq. (3.13.4) p 42]):

$$ \frac{\sigma(z+w)\sigma(z-w)}{\sigma(z)^2\sigma(w)^2} = \rho(w) - \rho(z) $$

$$ \frac{\sigma(2z)}{\sigma(z)^4} = -\rho'(z) $$

$$ \frac{\sigma(x+y+z)\sigma(x)\sigma(y)\sigma(z)}{\sigma(x+y)\sigma(y+z)\sigma(x+z)} = \frac{-1}{2} \begin{vmatrix} 1 & \wp(x) & \wp'(x) \\ 1 & \wp(y) & \wp'(y) \\ 1 & \wp(z) & \wp'(z) \end{vmatrix} \frac{1}{(\wp(x) - \wp(y))(\wp(y) - \wp(z))(\wp(z) - \wp(x))} $$

These allows to compute the cubical doubling, cubical differential additions and cubical three way additions on a Weierstrass model. There also exists Frobenius-Stickelberger relations for multiway relations (which allows to compute them directly rather than through three way additions), but we won't need them.

**Remark 4.35** (The normalised neutral point). Our cubical analytic neutral point $\tilde{0}$ is the one corresponding to $0 \in \mathbb{C}$. The normalisation is given by $(\sigma/(\wp/\wp'))(0) = -2$ by [OKUO11, Lemma 1]. Indeed $(\wp\sigma^2)(0) = 1$ and $(\wp'\sigma^3)(0) = -2$. But the Weierstrass equation is $\wp'^2 = \wp^3 - g_2\wp - g_3$, so for a Weierstrass equation of the form $y^2 = x^3 + ax + b$, via $x = \wp, y = \wp'/2$, $\tilde{0}$ corresponds to the normalisation $(\sigma/(-x/y))(0) = 1$, which is the opposite of the normalisation we use in Remark 4.24 (alternatively: $\sigma = -Z_1$).

### 4.9.3. Algebraic elliptic cubical points.
Our analytic cubical laws translate in the following algebraic law (taking into account the change of variable from Remark 4.35 between analytic Weierstrass coordinates and algebraic Weierstrass coordinates): given $\tilde{P} = (x(P), y(P), Z_1(\tilde{P}))$, $\tilde{Q} = (x(Q), y(Q), Z_1(\tilde{Q}))$ and $\tilde{P} - \tilde{Q} = (x(P - Q), y(P - Q), Z_1(\tilde{P} - \tilde{Q}))$, we have the cubical doubling

$$ 2\tilde{P} = (x(2P), y(2P), Z_1(\tilde{P})^4 2y(P)), $$

and the cubical differential addition is $\tilde{P} + \tilde{Q} = (x(P + Q), y(P + Q), Z_1(\tilde{P} + \tilde{Q}))$ with

$$ Z(\tilde{P} + \tilde{Q})Z(\tilde{P} - \tilde{Q}) = Z_1(\tilde{P})^2 Z_1(\tilde{Q})^2 (x(Q) - x(P)). $$

And the three way addition reads:

$$ \frac{Z_1(P + \tilde{Q} + R)Z_1(\tilde{P})Z_1(\tilde{Q})Z_1(\tilde{R})}{Z_1(\tilde{P} + \tilde{Q})Z_1(\tilde{Q} + R)Z_1(\tilde{P} + R)} = \frac{\begin{vmatrix} 1 & x(P) & y(P) \\ 1 & x(Q) & y(Q) \\ 1 & x(R) & y(R) \end{vmatrix}}{(x(Q) - x(P))(x(R) - x(P))(x(R) - x(Q))} = \frac{l_{P,Q}(R)}{(x(R) - x(P))(x(R) - x(Q))'} $$

where $l_{P,Q}: y - ax - \beta$ is the equation of the line going through $P$ and $Q$.

Finally, for general cubical points in Riemann position $[\tilde{P}_1, \tilde{P}_2, \tilde{P}_3, \tilde{P}_4; \tilde{Q}_1, \tilde{Q}_2, \tilde{Q}_3, \tilde{Q}_4]$, if we let $P_1 = P'_1 + P'_2$, $P_2 = P'_1 - P'_2$, $P_3 = P'_3 + P'_4$, $P_4 = P'_3 - P'_4$, then $2P'_1 = P_1 + P_2$, $2P'_2 = P_1 - P_2$, $2P'_3 = P_3 + P_4$, $2P'_4 = P_3 - P_4$, and we have: $\tilde{Q}_1 = P'_3 - P'_2$, $\tilde{Q}_2 = P'_3 + P'_4$, $\tilde{Q}_3 = P'_1 - P'_4$, $\tilde{Q}_4 = P'_1 + P'_4$. We can then use the cubical differential additions above to obtain

$$ \frac{Z_1(\tilde{P}_1)Z_1(\tilde{P}_2)Z_1(\tilde{P}_3)Z_1(\tilde{P}_4)}{Z_1(\tilde{Q}_1)Z_1(\tilde{Q}_2)Z_1(\tilde{Q}_3)Z_1(\tilde{Q}_4)} = \frac{(x(P'_2) - x(P'_1))(x(P'_4) - x(P'_3))}{(x(P'_2) - x(P'_3))(x(P'_4) - x(P'_1))}. $$
For the reader who do not like analytic arguments, let us give a purely algebraic proof of the cubical arithmetic for elliptic curves in Weierstrass form. Recall from Remark 4.24 that for an elliptic curve in Weierstrass form $y^2 = x^3 + a_2x^2 + a_4x + a_6$, we take $\tilde{0}$ corresponding to the uniformiser $\pi_{0_E} = x/y$, i.e. $(Z_1/(x/y))(\tilde{0}) = 1$

By Example 4.3 and Section 4.5.2, for a differential addition, we have
$$ \frac{Z_1(P+\tilde{Q}+R)Z_1(P-\tilde{Q}+R)Z_1(\tilde{R})Z_1(\tilde{\tilde{R}})}{Z_1(-\tilde{Q}+R)Z_1(\tilde{Q}+R)Z_1(\tilde{P}+R)Z_1(\tilde{P}+\tilde{R})} = \gamma(R). $$
Here by Example 4.3, $\gamma = \frac{g_{\tilde{Q},-\tilde{Q}}(\cdot)}{g_{\tilde{Q},-\tilde{Q}}(P-\cdot)}$.

For a generic $\tilde{Q}$, we can take $g_{\tilde{Q},-\tilde{Q}} = x - x(\tilde{Q})$ and so $\gamma(R) = \frac{x(P-R)-x(\tilde{Q})}{x(R)-x(\tilde{Q})}$. We also have $Z_1(-\tilde{Q}) = -Z_1(\tilde{Q})$ and even $Z_1(-\tilde{\tilde{Q}}) = -Z_1(\tilde{\tilde{Q}})$.

Since our choice of base point evaluation to represent cubical functions is $\tilde{0}$, we have
$$ Z_1(P+\tilde{Q})Z_1(P-\tilde{Q}) = Z_1(\tilde{P})^2Z_1(\tilde{Q})Z_1(-\tilde{Q})Z_1(\tilde{\tilde{Q}})^{-2}\gamma(0_E). $$

And since $Z_1/(x/y)(\tilde{0}) = 1$, we have $(\gamma/Z_1^2)(\tilde{0}) = (\gamma y^2/x^2)(0_E)$. Hence, we have
$$ Z_1(P+\tilde{Q})Z_1(P-\tilde{Q}) = Z_1(\tilde{P})^2Z_1(\tilde{Q})^2 - (\gamma y^2/x^2)(0_E). $$

The numerator of $\gamma$ evaluated at $R = 0_E$ gives $x(P) - x(\tilde{Q})$. The denominator of $\gamma$ multiplied by $(x/y)^2$ is $x(R)^3/y(R)^2 - x(R)^2x(\tilde{Q})/y(R)^2$. The term $x(R)^3/y(R)^2$ is equal to 1 at $R = 0_E$ because of the curve equation, while the second term gives 0. So $(\gamma y^2/x^2)(0_E) = x(P) - x(\tilde{Q})$. In the end, we have $Z_1(P+\tilde{Q})Z_1(P-\tilde{Q}) = (x(\tilde{Q}) - x(P))Z_1(\tilde{P})^2Z_1(\tilde{\tilde{Q}})^2$, as expected.

We let the reader work out the case of doublings and three way additions.

### 4.9.4. Cubical arithmetic on Kummer lines.
We could use $Z_2 = Z_1^2$ to compute the cubical arithmetic with respect to the divisor $D = 2(0_E)$: we would represent $\tilde{P}$ by $(x(P), y(P), Z_1^2(P))$. This combines level 3 projective coordinates with one level 2 affine coordinates.

Instead, we will use level 2 affine coordinates only, at the cost of moving to the Kummer line. The projective coordinates with respect to $2(0_E)$ are $(X_2 : Z_2)$, and so the associated affine representation from Section 4.5 is $\tilde{P} = (X_2(\tilde{P}), Z_2(\tilde{P})) = (x(P)Z_2(\tilde{P}), Z_2(\tilde{P}))$. It is enough to recover $P$ up to a sign, hence we can interpret this latter representation as a cubical representation over the Kummer line.

The cubical arithmetic in this affine representation then becomes, for cubical doublings, $2\tilde{P} = (x(2P)Z_2(2\tilde{P}), Z_2(2\tilde{P}))$ with $Z_2(2\tilde{P}) = 4y(P)^2Z(\tilde{P})^4$, and we can use the Weierstrass equation to express $y(P)^2$ in terms of $x(P)$. And the cubical differential addition is
$$ P+\tilde{Q} = (x(P+\tilde{Q})Z_2(P+\tilde{Q}), Z_2(P+\tilde{Q})) \text{ with } Z_2(P+\tilde{Q})Z_2(P-\tilde{Q}) = Z_2(\tilde{P})^2Z_2(\tilde{\tilde{Q}})^2(x(\tilde{Q}) - x(P))^2. $$

The three way addition is more complex. First from $x(P_1), x(P_2), x(P_1 + \tilde{Q}), x(P_2 + \tilde{Q})$ we can use a compatible addition to recover $x(P_1 - P_2)$ or $x(P_1 + P_2 + \tilde{Q})$, and then use a projective differential addition to recover $x(P_1 + P_2)$. Now we can combine the level 1 cubical differential additions formulas for $Z_1(P_1 + \tilde{P}_2 + \tilde{Q})Z_1(P_1 - \tilde{P}_2 + \tilde{Q})$, $Z_1(P_1 + \tilde{P}_2 + \tilde{Q})Z_1(P_2 - \tilde{P}_1 + \tilde{Q})$, $Z_1(P_1 - \tilde{P}_2 + \tilde{Q})Z_1(P_2 - \tilde{P}_1 + \tilde{Q})$, and $Z_1(P_1 + P_2)Z_1(P_1 - P_2)$, to obtain:

$$ 
\begin{aligned}
Z_2(P_1 + \tilde{P}_2 + \tilde{Q}) &= \frac{Z_2(P_1 + \tilde{Q})Z_2(P_2 + \tilde{Q})Z_2(\tilde{P}_1)Z_2(\tilde{P}_2)}{Z_2(P_1 - \tilde{P}_2)Z_2(\tilde{Q})} \frac{(x(P_2) - x(P_1 + \tilde{Q}))(x(P_1) - x(P_2 + \tilde{Q}))}{(x(P_1 - \tilde{P}_2) - x(\tilde{Q}))} \\
&= \frac{(X_2(\tilde{P}_2)Z_2(P_1 + \tilde{Q}) - X_2(P_1 + \tilde{Q})Z_2(\tilde{P}_2))(X_2(\tilde{P}_1)Z_2(P_2 + \tilde{Q}) - X_2(P_2 + \tilde{Q})Z_2(\tilde{P}_1))}{(X_2(P_1 - \tilde{P}_2)Z_2(\tilde{Q}) - Z_2(P_1 - \tilde{P}_2)X_2(\tilde{Q}))}
\end{aligned}
$$

Since compatible additions give $x(P_1 + P_2 + \tilde{Q})$, knowing $Z_2(P_1 + \tilde{P}_2 + \tilde{Q})$ is enough to recover $X_2(P_1 + \tilde{P}_2 + \tilde{Q})$, hence $P_1 + \tilde{P}_2 + \tilde{Q}$. (We will see a more direct formula
for $X_2(P_1 + \widetilde{P}_2 + Q)$ in the Montgomery model in Section 5.2). Likewise, we can compute $P_1 + \widetilde{P}_2 = \text{DiffAdd}(\widetilde{P}_1, \widetilde{P}_2, P_1 - \widetilde{P}_2)$, then the biextension element corresponding to $[P_1 + \widetilde{P}_2, \widetilde{Q}; \tilde{0}, P_1 + \widetilde{P}_2 + Q]$ is equal to $g_{P_1, Q} *_1 g_{P_2, Q}$.

For the compatible additions, let $\kappa_{00} = X(P_1+P_2)X(P_1-P_2)$, $\kappa_{01} = X(P_1+P_2)Z(P_1-P_2) + X(P_1-P_2)Z(P_1+P_2)$, $\kappa_{11} = Z(P_1+P_2)Z(P_1-P_2)$, so that $(X(P_1+P_2) : Z(P_1+P_2)), (X(P_1-P_2) : Z(P_1-P_2))$ are solutions of the homogeneous system $P(X, Z) = \kappa_{11}X^2 - \kappa_{01}XZ + \kappa_{00}Z^2 = 0$. By symmetry, we can always write the $\kappa$ as biquadratic polynomials in $(X(P_1) : Z(P_1)), (X(P_2) : Z(P_2))$. Define similarly $\kappa'$ for $P'_1 = P_1 + Q, P'_2 = P_2+Q$ to obtain a polynomial $P'(X, Z)$. The point $(X(P_1-P_2) : Z(P_1-P_2))$ is a root of both $P(X, Z)$ and $P'(X, Z)$, and can be written as $(\kappa'_{01}\kappa_{00} - \kappa_{01}\kappa'_{00} : \kappa'_{11}\kappa_{00} - \kappa_{11}\kappa'_{00})$.

We can use the $\mathbb{G}_m$-action from Lemma 4.15 to get rid of divisions: the cubical points $[\lambda \cdot P_1 + P_2, \widetilde{Q}; \tilde{0}, \lambda \cdot P_1 + \widetilde{P}_2 + Q]$ and $[P_1 + P_2, \widetilde{Q}; \tilde{0}, P_1 + \widetilde{P}_2 + Q]$ represent the same biextension elements. This allows to replace divisions by multiplication.

We recall that a biextension point is then represented via $g_{P,Q} = [\widetilde{P}, \widetilde{Q}; \tilde{0}, P + Q]$. We could imagine a mix of the $Z_2 = Z_1^2$ representation and the affine lift representation: when computing $g_{P,Q}^{*1,\ell} = [\ell\widetilde{P}, \widetilde{Q}; \tilde{0}, \ell\widetilde{P} + Q]$, compute $(X_2(\ell\widetilde{P}), Z_2(\ell\widetilde{P}))$ via a cubical biextension ladder to represent $\ell\widetilde{P}$, and use $Z_2(\ell\widetilde{P} + Q) = Z_2(\ell\widetilde{P} + Q)$ to represent $\ell\widetilde{P} + Q$. Indeed, to recover $X_2(\ell\widetilde{P} + Q)$, given $Z_2(\ell\widetilde{P} + Q)$ (and provided it is not zero), we just need to compute $x_2(\ell\widetilde{P} + Q)$. Since we are doing a ladder, we have $(\ell - 1)P$ and $\ell P$, and we can recover $\ell P + Q$ via a compatible addition $\text{CompatibleAdd}(\ell P, Q; (\ell - 1)P, P + Q)$.

So $Z_2(\ell\widetilde{P} + Q)$ is enough to compute the cubical arithmetic. Can we find fast formulas using $Z_2$ only without recovering $X_2$? This is related to the question of denominator elimination in Section 6.1.

### 4.9.5. Elliptic nets.
Elliptic nets are another way to compute the cubical arithmetic of an elliptic curve (or abelian variety), associated to the principal polarisation $(0_E)$. In elliptic nets, the cubical point $\widetilde{P} = (P, \sigma(\widetilde{P}))$ is only represented through $\sigma(\widetilde{P})$. This is not enough to determine $\widetilde{P}$ let alone do cubical arithmetic.

The key insight of Stange is the following: while $\sigma(\widetilde{P})$ alone is not enough, the data of $\sigma(m\widetilde{P})$ for small values of $m$ is enough to recover any $\sigma(n\widetilde{P})$ through a recurrence formula. Likewise, the values of $\sigma(n\widetilde{P} + mQ)$ for small values of $n, m$ allow to compute all values $\sigma(\ell_1\widetilde{P} + \ell_2Q)$ through a more complicated recurrence formula.

The resulting recurrence formulas allows to build iteratively the so called elliptic nets. Since elliptic nets compute the cubical arithmetic, they can be used to compute pairings: compare [Sta08, Theorems 17.2.1, 17.2.2] with Theorem 4.19 and [OKU01, Theorems 4, 5] with Proposition 4.22. A somewhat annoying thing with the elliptic net representation is that it cannot represent $\widetilde{0}_E$ since $\sigma(0_E) = 0$, so when using it to compute pairings through the bixtension monodromy, it requires to compare $g_{P,Q}^{\ell+1}$ with $g_{P,Q}$ rather than $g_{P,Q}^{\ell}$ with the constant 1.

Elliptic nets extend to abelian varieties using analytic thetas instead of the Weierstrass $\sigma$ function, see [Tra14, Chapitre 3] (which also gives formula for analytic thetas of higher level than $n=2$). Although the result is stated there for Jacobians of hyperelliptic curves, the proof is valid for any abelian variety. The idea is that the differential addition formula deduced from Riemann relations give a bilinear relation between products of the form $\theta_i(z_1+z_2)\theta_j(z_1-z_2)$ and products of the form $\theta_i(z_1)\theta_i(z_2)$ and $\theta_i(z_2)\theta_j(z_2)$. The associated matrix to these bilinear relations is not of full rank and taking its determinant give the recurrence relation sought for.
We now compare elliptic nets with our approach. First the strategy of the proof is slightly different: Stange uses the analytic $\sigma$ function to derive algebraic formulas for elliptic nets over $\mathbb{C}$, which she uses to give algebraic formulas for elliptic nets over any field. Here we use the fact that we already know that the cubical arithmetic is algebraic to give an algebraic interpretation of the analytic function $\sigma$ as defining a cubical point, so that we can define an algebraic $\sigma$ over any field. We could then recover the algebraic elliptic net recurrence from this algebraic $\sigma$.

But the main difference is in terms of our choice of representation of cubical point: representing them only by $\sigma(\tilde{P})$ loses a lot of information and the elliptic net representation needs to work with many values. In Section 4.9, we give cubical formulas for the representation $(P, \sigma(\tilde{P}))$ where we keep track of the underlying point on the elliptic curve together with $\sigma(\tilde{P})$. In other words: the cubical points are the intrinsic object, and since $\sigma(\tilde{P})$ alone can not recover $\tilde{P}$, it makes working with the elliptic nets $\sigma(\sum n_i \tilde{P}_i)$ directly harder than working with the cubical points $\sum n_i \tilde{P}_i$ themselves.

And we will see in Section 5 that the sweet spot for efficient cubical arithmetic seems to be to use the affine representation associated to the divisor $2(0_E)$ and the global sections $X, Z$, i.e. represent $\tilde{P}$ via the level two affine coordinates $X(\tilde{P}), Z(\tilde{P})$ rather than by only one affine coordinate $\sigma(\tilde{P})$. Via this representation, as explained in Section 4.5 we will be able to leverage the efficient Kummer arithmetic to an efficient cubical arithmetic.

The only drawback is that the polarisation associated to $2(0_E)$ is not principal anymore, so as explained in Remark 2.11 we will need to use the full power of Theorem 2.9 rather than Corollary 2.5 to correctly handle the case $\ell$ even.

## 5. PAIRINGS ON KUMMER LINES

We specialize in Section 5.1 our whole framework of cubical arithmetic from Section 4 to the case of Kummer lines. Then we specialize further in Section 5.2 to the case of the Montgomery model of Kummer lines, which allows us to prove Theorem 1.1.

5.1. **Algorithms for a Kummer line.** We will represent our cubical points through the affine lift representation associated to the sections $X, Z \in \Gamma(2(0_E))$ of the divisor $D = 2(0_E)$. We will always take $\tilde{0} = (1, 0)$ as our choice of affine lift of $0_E = (1:0)$. Indeed, we have seen in Remark 4.34 that normalizing $\tilde{0}$ in level 2 through $(Z_2/(x/y)^2)(\tilde{0}) = 1$, makes $X_2(\tilde{0}) = 1$, so $\tilde{0} = (1, 0)$ in level 2 corresponds to $Z_2(\tilde{0}) = 1$ being normalised with respect to $(x/y)^2$.

The algorithm to compute the non reduced Tate pairing is given in Algorithm 5.1. To compute the Weil pairing, it suffices to apply Algorithm 5.1 again with the arguments $P, Q$ reversed, but making the same choices for $\tilde{P}, \tilde{Q}, P + Q$ (or if making different choices for whatever reason, ensuring the associated biextension element $g_{Q,P}$ is equal to $\ell(g_{P,Q})$, to obtain $\lambda_Q$. The Weil pairing is $e_{W,2(0_E),\ell}(P, Q) = \lambda_P/\lambda_Q$.

This is the (non reduced) Tate and Weil pairings associated to $2(0_E)$, so the square of the usual non reduced Tate and Weil pairing. When $\ell$ is odd, and if the characteristic $p > 0$, to recover the usual reduced Tate pairing it suffices to adjust the final exponentiation from $(q-1)/\ell$ to $(q-1)/(2\ell)$. In any case, having the squares of the usual pairings on $\mu_\ell$ is not a problem in practice if $\ell$ is odd (one could always take a square root if needed). This becomes a problem if $\ell$ is even, in which case we use the strategy of Remark 2.11 to compute the standard (non reduced) Tate and Weil pairing associated to $(0_E)$. In this case, $\ell = 2m$, and rather than computing $\ell \tilde{P}, \ell \tilde{P} + Q$, we compute $\tilde{m} \tilde{P}, \tilde{m} \tilde{P} + Q$.

By assumption, $T = mP$ is a point of two torsion, and we can now use the action of an element of the theta group $g_T \in G(2(0_E))$ above $T$ as explained in Section 4.5.
**Input:** $(X(P) : Z(P))$ a rational point of $\ell$-torsion, $(X(Q) : Z(Q))$ a rational point, $(X(P+Q) : Z(P+Q))$

**Output:** The non reduced Tate pairing $e_{T,2(0_E),\ell}(P,Q)$

→ Take arbitrary affine lifts $\tilde{P} = (X(\tilde{P}), Z(\tilde{P}))$, $\tilde{Q} = (X(\tilde{Q}), Z(\tilde{Q}))$, $\overline{P+Q} = (X(\overline{P+Q}), Z(\overline{P+Q}))$.

→ Compute $\ell\overline{P}, \ell\overline{P+Q}$ using either Algorithms 4.2 and 4.7.

→ Write $\overline{\ell P} = \lambda_{0,P}\tilde{0}, \ell\overline{P+Q} = \lambda_{1,P}\tilde{Q}$.

→ Return $\lambda_P = \lambda_{1,P}/\lambda_{0,P}$.

<caption>ALGORITHM 5.1. The non reduced Tate pairing associated to $2(0_E)$ on a Kummer line</caption>

Concretely, the translation by $T$ is given by a projective matrix, i.e. an action of the form $(X, Z) \mapsto (aX + bZ : cX + dZ)$. Take any choice of affine lift $(X, Z) \mapsto (aX + bZ, cX + dZ)$ of the translation by $T$ on $(X, Z)$ coordinates. In other words, we lift the projective $2 \times 2$ matrix of translation to a standard matrix; this is the same as making a choice of theta group element $g_T \in G(2(0_E))$ above $T$ and looking at the action of $g_T$ on $(X, Z)$. The choice of an affine point $\tilde{T}$ is a way to encode this choice: $g_T$ should then be the unique element such that $g_T \cdot \tilde{0} = \tilde{T}$. We call the action of a $g_T$ corresponding to $\tilde{T}$ on a cubical point $\tilde{P}$ the translation of $\tilde{P}$ by $\tilde{T}$ and denote it $\overline{\tilde{P} + \tilde{T}} = \text{Translate}(\tilde{P}, \tilde{T})$. Then we can compute $\overline{\ell \tilde{P}} = \overline{m \tilde{P} + \tilde{T}}, \ell \overline{\tilde{P} + Q} = m \overline{\tilde{P} + Q} + \tilde{T}$.

Since $T = mP$, we could use in practice the choice of $\tilde{T} = \overline{mP}$ which just computed; but we remark that a different choice of $\tilde{T}$ would still get the correct value as long as we use the same choice of $\tilde{T}$ for both the translation action to compute $\overline{\ell \tilde{P}}$ and $\ell \overline{\tilde{P} + Q}$, since changing $\tilde{T}$ to $\lambda \tilde{T}$ would change $\overline{\ell \tilde{P}}$ to $\lambda \overline{\ell \tilde{P}}$ and $\ell \overline{\tilde{P} + Q}$ to $\lambda \ell \overline{\tilde{P} + Q}$, so they would still induce the same biextension element by Lemma 4.15.

We will also denote by Translate$(\tilde{P}, T)$ an algorithm which makes a choice of $\tilde{T}$ depending on $T$ and then call Translate$(\tilde{P}, \tilde{T})$; as long as the choice of $\tilde{T}$ is the same for every call with the same $T$, the computation will be correct. This can save us some arithmetic operations, for instance on a Montgomery model if $\overline{m\tilde{P}} = (0, a)$, we prefer to take $(\tilde{0} : 1) = (0, 1)$ rather than $(0, a)$. Indeed in the former case, the translation action by $(0, 1)$ is simply $(X, Z) \mapsto (Z, X)$, while the translation action of $(0, a)$ would be $(X, Z) \mapsto (aZ, aX)$, so we save multiplications.

The algorithm for the non reduced Tate pairing for even $\ell$ is given in Algorithm 5.2, and like for Algorithm 5.1 it gives an algorithm for the Weil pairing by calling it a second time with $(P, Q)$ swapped and computing $e_{W,(0_E),\ell}(P, Q) = \lambda_P/\lambda_Q$.

By Remark 4.27, it is important in Algorithm 5.1 that we normalize our cubical points via $Z = 1$. Indeed, the divisor of zeroes of $Z$ is $2(0_E)$, so it satisfies the condition of Algorithm 5.1, and normalising all points with $Z = 1$ is a convenient way to be sure that the associated biextension element $g_{P,Q}$ is a rational tensor square. (We do have $Z(\tilde{0}) = 0$, but $\tilde{0}$ is normalised for $Z$ with respect to $(x/y)^2$ which is a square too.)

We could allow an arbitrary normalisation of $\tilde{P}, \tilde{Q}, \overline{P+Q}$, but we would then need to adjust $\lambda_P$ by $\left( \frac{Z(\overline{P+Q})}{Z(\overline{P})Z(\overline{Q})} \right)^{\ell}$.
**Input:** $(X(P) : Z(P)) \text{ a rational point of } \ell\text{-torsion}, (X(Q) : Z(Q)) \text{ a rational point},$ 
$(X(P+Q) : Z(P+Q))$

**Output:** The non reduced Tate pairing $e_{T,(0_E),\ell}(P, Q)$

-> Take affine lifts $\tilde{P} = (X(\tilde{P}), Z(\tilde{P})) = 1$, $\tilde{Q} = (X(\tilde{Q}), Z(\tilde{Q})) = 1$, $\overline{P+Q} = (X(\overline{P+Q}), Z(\overline{P+Q})) = 1$.

-> Compute $\overline{m\tilde{P}}, m\tilde{P} + \tilde{Q}$ using either Algorithms 4.2 and 4.7.

-> Let $T = mP$, and compute $\overline{\ell P} = \text{Translate}(\overline{m\tilde{P}}, T)$ and $\overline{\ell P + Q} = \text{Translate}(m\tilde{P} + Q), T)$.

-> Write $\overline{\ell P} = \lambda_{0,P}\tilde{0}, \overline{\ell P + Q} = \lambda_{1,P}\tilde{Q}$.

-> Return $\lambda_P = \lambda_{1,P} / \lambda_{0,P}$.

ALGORITHM 5.2. The non reduced Tate pairing associated to (0_E) on a Kummer line for l = 2m even

**Input:** $(X(P) : Z(P)) \in \mathbb{G}_2, (X(Q) : Z(Q)) \in \mathbb{G}_1, (X(P+Q) : Z(P+Q)), \lambda \equiv q \mod \ell$

**Output:** The non reduced Ate pairing $a_{\lambda,2(0_E),\ell}(P, Q)$

-> Take arbitrary affine lifts $\tilde{P} = (X(\tilde{P}), Z(\tilde{P}))$, $\tilde{Q} = (X(\tilde{Q}), Z(\tilde{Q}))$, $\overline{P+Q} = (X(\overline{P+Q}), Z(\overline{P+Q}))$.

-> Compute $\overline{m\tilde{P}}, m\tilde{P} + \tilde{Q}$ using either Algorithms 4.2 and 4.7.

-> Write $\overline{\lambda P} = \lambda_{0,P}\pi_q(\tilde{P}), \overline{\lambda P + Q} = \lambda_{1,P}\pi_q(\tilde{P} + Q)$.

-> Return $\lambda_P = \lambda_{1,P} / \lambda_{0,P}$.

ALGORITHM 5.3. The non reduced Ate pairing associated to 2(0<sub>E</sub>) on a Kummer line

We can also use Proposition 4.22 to compute the (non reduced) Ate pairing (with respect to 2(0<sub>E</sub>), so the square of the usual Ate pairing). We reuse the notations of Section 3.4 and Proposition 4.22, and let $m \equiv \ell \mod q$. We recall that $\pi_q((X(\tilde{P}), Z(\tilde{P}))) = (\pi_q(X(\tilde{P})), \pi_q(Z(\tilde{P})))$.

A summary is given by:

**Theorem 5.1.** Let $E/\mathbb{F}_q$ be an elliptic curve, and $X = X_{2(0_E)}$ the biextension associated to the divisor $2(0_E)$.

When $\ell$ is odd, to compute the square of the non reduced Tate pairing (resp. of the Weil pairing), we need to compute one (resp. two) biextension exponentiation by $\ell$.

When $\ell$ is even, to compute the usual non reduced Tate pairing (resp. of the Weil pairing), we need to compute one (resp. two) biextension exponentiation by $\ell/2$, followed by an affine translation by a point of two torsion.

To compute the square of the usual $m$-Ate pairing (with $m \equiv q \mod \ell$), we need to compute one biextension exponentiation by $m$.
Furthermore, using the cubical representation of biextension elements, and representing cubical points $\tilde{P}$ by the affine lifts $\tilde{P} = (X(\tilde{P}), Z(\tilde{P}))$, where $X, Z$ are sections of $2(0_E)$, a biextension exponentiation then costs:

* Using the cubical biextension ladder from Algorithm 4.2, one affine doubling and two affine differential additions by bits. (One could replace an affine differential addition by a three way addition instead.)

As special cases:

- When $l = 2^m$ or for self pairing $e(P, P)$, one affine doubling and one affine differential addition by bits

- For batch Tate or Ate pairings computations $e(P, Q_i)$ with the same base point $P$, after the first biextension exponentiation the follow up biextension exponentiation only cost one differential addition (alternatively, one three way addition) by bits

* Using the cubical biextension double and add from Algorithm 4.1, one affine doubling and one differential addition (or alternatively one three way addition) for each doubling, and one compatible addition and one three way addition for each addition.

### 5.2. The Montgomery model.
As explained in Section 5.1, to compute pairings using the formulas from Theorem 4.19 and Proposition 4.22, we need to compute cubical biextension exponentiations using either Algorithm 4.2 or Algorithm 4.7.

#### 5.2.1. Cubical ladder on the Montgomery model.
For the first one, we just need to explain how to compute affine doublings and affine differential additions in $X, Z$ coordinates. For the second one, we also need to give algorithms for compatible additions and three way additions in $X, Z$ coordinates.

In this section, we thus concentrate on the Montgomery model of Kummer lines. We use Section 4.9 to derive our cubical arithmetic. We start with a Montgomery curve $By^2 = x^3 + Ax^2 + 1$. $x$-only additions have a particularly nice form on the Montgomery model:

$$x(2P) = \frac{(x(P)^2 - 1)^2}{4x(P)(x(P)^2 + Ax(P) + 1)} = \frac{(X(P)^2 - Z(P)^2)^2}{4X(P)Z(P)(X(P)^2 + AX(P)Z(P) + Z(P)^2)},$$

and

$$x(P+Q)x(P-Q) = \left( \frac{x(P)x(Q) - 1}{x(P) - x(Q)} \right)^2 = \left( \frac{X(P)X(Q) - Z(P)Z(Q)}{X(P)Z(Q) - Z(P)X(Q)} \right)^2.$$

On the other hand, by Section 4.9.4, we have

$$Z(2\tilde{P}) = 4X(\tilde{P})Z(\tilde{P})(X(\tilde{P})^2 + AX(\tilde{P})Z(\tilde{P}) + Z(\tilde{P})^2),$$

and

$$Z(P+Q)Z(P-Q) = (X(Q)Z(\tilde{P}) - X(\tilde{P})Z(Q))^2.$$

In particular, combining both equations, we see that the natural way to write $x(2P)$, $x(P+Q)x(P-Q)$ as rational functions in terms of the projective coordinates $X(P), Z(P), X(Q), Z(Q)$ already gives us the correct cubical arithmetic by taking the numerator for $X(2\tilde{P}), X(P+Q)X(P-Q)$ and the denominator for $Z(2\tilde{P}), Z(P+Q)Z(P-Q)$ respectively! We remark that we can also directly recover the formulas for $X(2\tilde{P})$ and $X(P+Q)X(P-Q)$ from the ones for $Z$, by using that if $\tilde{T} = (1, 0)$, $X(\tilde{P}) = Z(P+T)$. We obtain Algorithms 5.4 and 5.5 for the cubical arithmetic.

Some comments are in order: in these algorithms we take $\tilde{0} = (1, 0)$. Taking a different affine lift of $\tilde{0}$, we would need to adjust the doubling and differential addition formulas accordingly. Also the Montgomery coefficients $\mathcal{A}$ or $\frac{\mathcal{A}+2}{4}$ are often represented in projective
Input: $\tilde{P} = (X(P), Z(P))$

Output: $\tilde{2P} = (X(2P), Z(2P)) = \text{Double}(\tilde{P})$

$$ \rightarrow a = (X(P) + Z(P))^2 $$

$$ \rightarrow b = (X(P) - Z(P))^2 $$

$$ \rightarrow c = a - b $$

$$ \rightarrow X(2P) = ab $$

$$ \rightarrow Z(2P) = c(b + \frac{\lambda+2}{4}c) $$

ALGORITHM 5.4. Affine cubical doubling in the Montgomery model

Input: $\tilde{P} = (X(P), Z(P)), \tilde{Q} = (X(Q), Z(Q)), P - \tilde{Q} = (X(P - Q), Z(P - Q))$

Output: $P + \tilde{Q} = (X(P + Q), Z(P + Q)) = \text{DiffAdd}(\tilde{P}, \tilde{Q}, P - \tilde{Q})$

$$ \rightarrow u = (X(P) + Z(P))(X(Q) - Z(Q)) $$

$$ \rightarrow v = (X(P) - Z(P))(X(Q) + Z(Q)) $$

$$ \rightarrow 4X(P + Q)X(P - Q) = (u + v)^2 $$

$$ \rightarrow 4Z(P + Q)Z(P - Q) = (u - v)^2 $$

ALGORITHM 5.5. Affine cubical differential addition in the Montgomery model

coordinates $\frac{\lambda+2}{4} = (A_{24} : C_{24})$. Here I gave the cubical formulas where $C_{24}$ is normalised to 1. (Alternatively, using the notations from Algorithm 5.4, a $C_{24}$ not normalised to 1 corresponds to doubling formulas given by $X(2P) = C_{24}ab, Z(2P) = c(C_{24}b + A_{24}c)$, and associated to a different normalisation of $\tilde{0}$, but then we would need to update the cubical differential additions for this different normalisation.) However, as explained in Remark 4.23, for pairings we only need the biextension arithmetic, which is more flexible than the cubical arithmetic. In particular, the extra factor 4 in the affine differential addition formulas do not matter for pairings (it does for other applications though): in the case of the Weil pairing they are compensated in the quotient $\lambda_P/\lambda_Q$, and in the case of the Tate pairing it is killed by the final exponentiation, as long as we are not over the base field $\mathbb{F}_p$.

To compute pairings for $\ell$ even, we also need to make a choice of affine lift of translation by points of 2-torsion (i.e, fix once and for all an element $(T, g_T) \in G(2(0_E))$ represented via its action on the $(X, Z)$ coordinates). This is done in Algorithm 5.6. We remark also that for Montgomery Kummer lines, we could apply Algorithm 5.2 by normalising our points with respect to $X$ rather than to $Z$, because the divisor of zeroes of $X$ is $2(T)$, with $T = (0 : 1)$, so by Remark 4.27, $X$ is also a suitable coordinate for normalisation.

**Example 5.2.** Let $P$ be a point of 2-torsion. We want to compute the Tate pairing $e_{T,2}(P, Q)$ with respect to $(0_E)$ rather than with respect to $2(0_E)$ (since the latter would be trivial). By Remark 2.11, we can use either the $X$ or $Z$ coordinate, since the divisor of $Z$ is $2(0_E)$ and the divisor of $X$ is $2(T)$, with $T = (0 : 1)$, so they are both twice a rational divisor.

Let $M_P$ be the matrix of translation by $P$ in $(X, Z)$ coordinates (formally this is the matrix of some theta group element $g_P \in G(2(0_E))$ above $P$ action on $X, Z$). There are two approaches to this pairing computation:
**Input:** $\tilde{P} = (X(P), Z(P))$ and $T = (X(T) : Z(T))$ a point of 2-torsion

**Output:** $\tilde{P}+T = (X(P+T), Z(P+T)) = \text{Translate}(\tilde{P}, T)$

→ If $T = 0_E = (1:0)$, return $(X(P), Z(P))$

→ If $T = (0:1)$, return $(Z(P), X(P))$

→ Else $T = (X(T) : Z(T))$, we fix once and for all $\tilde{T} = (X(T), Z(T))$ and return $(X(T)X(P) - Z(T)Z(P), Z(T)X(P) - X(T)Z(P))$

<caption>ALGORITHM 5.6. Affine cubical translation by a point of 2-torsion in the Montgomery model</caption>

* We start with $\tilde{P}, \tilde{Q}, \tilde{P}+\tilde{Q}$ normalised to have $Z(\tilde{P}) = Z(\tilde{Q}) = Z(\tilde{P}+\tilde{Q}) = 1$ (alternatively normalized by their X coordinates). We compute the monodromy actions $M_P \cdot \tilde{P}+\tilde{Q} = \lambda_{P,1}\tilde{Q}$, and $M_P \cdot \tilde{P} = \lambda_{P,0}\tilde{0}$, and return the Tate pairing as $e_{T,2}(P, Q) = \lambda_{P,1}/\lambda_{P,0}$.

* We start with $\tilde{Q}$ arbitrary, but we use $\tilde{P} = M_P \cdot \tilde{0}, \tilde{P}+\tilde{Q} = M_P \cdot \tilde{Q}$. Since P is of two torsion, $M_P^2$ is a constant times the identity matrix: $M_P^2 = \lambda_P \text{Id}$. We have $M_P \cdot \tilde{P}+\tilde{Q} = \lambda \tilde{Q}, M_P \cdot \tilde{P} = \lambda \tilde{0}$, so the monodromy $\lambda/\lambda = 1$ is trivial. However, since our points $\tilde{P}, \tilde{Q}, \tilde{P}+\tilde{Q}$ are no longer normalized, we need to keep track of their coordinates to get the correct Tate pairing: $e_{T,2}(P, Q) = \frac{Z(\tilde{P})Z(\tilde{Q})}{Z(\tilde{P}+\tilde{Q})}$.

It is not hard to see that both methods give the same result; this is a special case of the more general formula from Remark 2.11 (which handle the case where the pairing is recovered from both a non trivial monodromy action and a non trivial normalisation).

We stress that these two methods work for $Z$ (or $X$), because their zero divisors are twice a divisor; it would not work for a random linear combination $U = u_0X + u_1Z$.

In practice, applying the second method, if $\tilde{P} = (X(P), Z(P)), \tilde{Q} = (X(Q), Z(Q))$, then $\tilde{P}+\tilde{Q} = (X(P)X(Q) - Z(P)Z(Q), Z(P)X(Q) - X(P)Z(Q))$. We obtain $e_{T,2}(P, Q) = \frac{Z(\tilde{P})Z(\tilde{Q})}{Z(\tilde{P}+\tilde{Q})} = \frac{Z(P)Z(Q)}{Z(P)X(Q)-X(P)Z(Q)} = \frac{1}{x(Q)-x(P)}$.

This is indeed the correct formula for the non reduced Tate pairing (recall that we take the opposite sign convention than the usual one): the normalised function $\mathbf{f}_{2,P} = x - x(P)$, so its evaluation at $Q$ is $x(Q) - x(P)$.

Using the X coordinate instead, we have $e_{T,2}(P, Q) = \frac{X(P)X(\tilde{Q})}{X(\tilde{P}+\tilde{Q})X(\tilde{0})} = \frac{X(P)X(Q)}{X(P)X(Q)-Z(P)Z(Q)} = \frac{1}{1-\frac{1}{x(P)x(\tilde{Q})}}$. And indeed, we have $\mathbf{f}_{2(P+T)-2(T)}((Q)-(0)) = \mathbf{f}_{2(P)-2(0)}((Q+T)-(T)) = \frac{x(Q+T)-x(P)}{x(T)-x(P)} = \frac{1/x(Q)-x(P)}{-x(P)} = 1-\frac{1}{x(P)x(Q)}$.

**Remark 5.3 (Complete arithmetic laws).** We remark that Algorithm 5.4 is always well defined, even for doubling $\tilde{0}$, and Algorithm 5.5 is also always well defined, except in the case where $\tilde{P}-\tilde{Q}$ is above a point $P-Q = 0_E = (0:1)$ or $P-Q = T = (1:0)$. In the first case, we have $P=Q$, and so we can first compute $2\tilde{Q} = \text{Double}(\tilde{Q})$ and then adjust the result depending on the scalar thus that $\tilde{P} = \lambda_P \cdot \tilde{Q}, \tilde{P}-\tilde{Q} = \lambda_0\tilde{0}$: $\text{DiffAdd}(\tilde{P}, \tilde{Q}, \tilde{P}-\tilde{Q}) = \frac{\lambda_P^2}{\lambda_0} \text{Double}(\tilde{Q})$. In the second case, we first lift the translation action by $T$ by taking $\tilde{T} = (1,0) = g_T \cdot \tilde{0}$, with $g_T \in G(2(0_E))$ satisfies $g_T \cdot (X,Z) = (Z,X)$. Now by Lemma 4.12,
$\text{DiffAdd}(g_T \cdot \tilde{Q}, \tilde{Q}, \tilde{T}) = g_T \cdot \text{Double}(\tilde{Q})$. We can thus compute $\text{Double}(\tilde{Q})$, apply $g_T$ on it (i.e. switch the X and Z coordinates). Then we need to adjust by the scalars such that
$$ \tilde{P} = \lambda_P \cdot g_T \cdot \tilde{Q}, P - \tilde{Q} = \lambda_0 \tilde{T}. \text{ Then } \text{DiffAdd}(\tilde{P}, \tilde{Q}, P - \tilde{Q}) = \frac{\lambda_P^2}{\lambda_0} g_T \cdot \text{Double}(\tilde{Q}). $$

This means that, adapting our formulas for these particular case, we have complete cubical arithmetic laws on Montgomery Kummer lines. In practice, for a pairing computation $e(P, Q)$, we can always use Algorithms 5.4 and 5.5, except when $P, Q$ or $P + Q$ is equal to the two torsion point $T = (0 : 1)$. We might as well treat these cases globally, or use Corollary 4.13.

Still, it is an interesting exercise to unravel the cubical formulas for these special cases, and check we recover the correct pairings. Let $T = (1 : 0)$, we fix $\tilde{T} = (1, 0)$, since $\tilde{0} = (0, 1)$ this means that the corresponding theta group element $g_T \in G(2(0_E))$ is the one that acts on $(X, Z)$ by $(X, Z) \mapsto (Z, X)$.

If $P + Q = T$, then $Q = T - P$, so by bilinearity $e(P, Q) = e(P, T)e(P, -P)$, so we only need to treat the case $P = T$ and $Q = T$. Furthermore, these cases are only interesting when $\ell = 2m$ is even, otherwise the pairings are trivial.

Let us first work out the case $P = T$. We can take $\tilde{P} = \tilde{T} = (1, 0)$, fix $\tilde{Q} = (X(Q), Z(Q))$ arbitrarily, and take $P + \tilde{Q} = \tilde{T} \cdot \tilde{Q} = (Z(Q), X(Q))$. Then $m\tilde{P} = (1, 0)$ if $m$ is odd, and $(0, 1)$ if $m$ is even. But in the odd case, for the monodromy we need to act by $-m\tilde{P} = (1, 0)$, so $-m\tilde{P} \cdot m\tilde{P} = \tilde{0}$. By Corollary 4.13, $m\tilde{P} + \tilde{Q} = \tilde{Q}$ if $m$ is even, and is equal to $P + \tilde{Q}$ if $m$ is odd. In both cases $-m\tilde{P} \cdot m\tilde{P} + \tilde{Q}$ is equal to $\tilde{Q}$. So the monodromy is trivial in this case. Since the divisor is $2(0_E)$, this monodromy gives the square of the Tate pairing $e_{T,2}(P, Q)$, which is indeed trivial since $P$ is of 2-torsion. By Remark 2.11, if we want to recover $e_{T,2}(P, Q)$, we need to correct the monodromy information by the coordinates $\left( \frac{Z(P+\tilde{Q})}{Z(\tilde{P})Z(\tilde{Q})} \right)^m = (X(Q)/Z(Q))^m = x(Q)^m$. We indeed have $\mathbf{f}_{2m,P}((Q) - (0)) = \mathbf{f}_{2,P}((Q) - (0))^m = (x(Q) - x(P))^m = x(Q)^m$.

Now for the case $Q = T$; we take $\tilde{Q} = (1, 0)$, and $P + \tilde{Q} = \tilde{Q} \cdot \tilde{P} = (Z(P), X(P))$. We compute $m\tilde{P}$ and $m\tilde{P} + \tilde{Q}$, then we act by $-m\tilde{P}$ which we can take to be equal to $-m\tilde{P}$. By Corollary 4.13, $m\tilde{P} + \tilde{Q} = \tilde{Q} \cdot \tilde{m}\tilde{P}$. If $g_{\tilde{Q}}, g_{-m\tilde{P}}$ are the theta group elements encoded by $\tilde{Q}, -m\tilde{P}$, it follows that the monodromy information we get is precisely the commutator pairing between $g_{\tilde{Q}}, g_{-m\tilde{P}}$, hence by Example 2.1 is equal to $e_{W,2}(-m\tilde{P}, \tilde{Q})$. In particular, for the Weil pairing $e_{W,2m}(P, \tilde{Q}) = e_{W,2}(m\tilde{P}, \tilde{Q})$, since the monodromy information on the preceding paragraph was trivial, it follows that the monodromy we compute with respect of both $P, \tilde{Q}$ do give the Weil pairing. For the Tate pairing however, the monodromy is not enough, and we need to correct by the cubical function $\left( \frac{Z(P+\tilde{Q})}{Z(\tilde{P})Z(\tilde{Q})} \right)^m = (X(P)/Z(P))^m = x(P)^m$. The end result we obtain for the Tate pairing via cubical arithmetic is $e_{T,2m}(P, \tilde{Q}) = x(P)^m e_{W,2}(m\tilde{P}, \tilde{Q})$. Comparing with the usual formula, we indeed have $e_{T,2m}(P, \tilde{Q}) = \mathbf{f}_{2m,P}((\tilde{Q}) - (0)) = e_{W,2m}(P, \tilde{Q})\mathbf{f}_{2m,\tilde{Q}}((\tilde{P}) - (0)) = e_{W,2}(m\tilde{P}, \tilde{Q})x(P)^m$.

We obtain the following complexity, which complements Theorem 5.1, for pairings in the Montgomery model of Kummer lines:

**Theorem 5.4.** Given the Kummer line Montgomery coordinates $(X(P) : Z(P))$ of $P, Q, P+Q$, first do a batch inversion to compute $x(P) = X(P)/Z(P)$, $x(Q), x(P+Q)$ and their inverses, and also compute $\frac{\lambda+2}{4}$. Take for cubical points $\tilde{0} = (1, 0)$, $\tilde{P} = (x(P), 1)$, $\tilde{Q} = (x(Q), 1)$, $P + \tilde{Q} = (x(P+Q), 1)$.

Then during the cubical biextension ladder, each affine doubling costs $2M + 2S + 1m_0$ and each affine differential addition costs $3M + 2S$.
**Input:** $P_1 = (X(P_1) : Z(P_1)), P_2 = (X(P_2) : Z(P_2)), P_1 + Q = (X(P_1 + Q) : Z(P_1 + Q)), P_2 + Q = (X(P_2 + Q) : Z(P_2 + Q))$
**Output:** $P_1 + P_2 = (X(P_1 + P_2) : Z(P_1 + P_2))$

$$ \rightarrow \kappa_{00} = (X(P_1)X(P_2) - Z(P_1)Z(P_2))^2 $$

$$ \rightarrow \kappa_{11} = (X(P_1)Z(P_2) - X(P_2)Z(P_1))^2 $$

$$ \rightarrow \kappa_{01} = 2((X(P_1)X(P_2) + Z(P_1)Z(P_2))(X(P_1)Z(P_2) + X(P_2)Z(P_1))) + 2AX(P_1)X(P_2)Z(P_1)Z(P_2)) $$

$$ \rightarrow \kappa'_{00} = (X(P_1 + Q)X(P_2 + Q) - Z(P_1 + Q)Z(P_2 + Q))^2 $$

$$ \rightarrow \kappa'_{11} = (X(P_1 + Q)Z(P_2 + Q) - X(P_2 + Q)Z(P_1 + Q))^2 $$

$$ \rightarrow \kappa'_{01} = 2((X(P_1 + Q)X(P_2 + Q) + Z(P_1 + Q)Z(P_2 + Q))(X(P_1 = Q)Z(P_2 + Q) + X(P_2 + Q)Z(P_1 + Q)) + 2AX(P_1 + Q)X(P_2 + Q)Z(P_1 + Q)Z(P_2 + Q)) $$

$$ \rightarrow X(P_1 - P_2) = \kappa'_{01}\kappa_{00} - \kappa_{01}\kappa'_{00} $$

$$ \rightarrow Z(P_1 - P_2) = \kappa'_{11}\kappa_{00} - \kappa_{11}\kappa'_{00} $$

$$ \rightarrow X(P_1 + P_2) = \frac{\kappa_{00}}{X(P_1 - P_2)} $$

$$ \rightarrow Z(P_1 + P_2) = \frac{\kappa_{11}}{Z(P_1 - P_2)} $$

<caption>ALGORITHM 5.7. Compatible additions in the Montgomery model</caption>

Theorem 1.1 follows from Theorems 5.1 and 5.4.

5.2.2. Double and add biextension algorithm for the Montgomery model. For the double and add algorithm, we need to give the formulas for the compatible addition. Using the notations of Section 4.9.4, we have

$$ \kappa_{00} = (X(P_1)X(P_2) - Z(P_1)Z(P_2))^2, $$

$$ \kappa_{11} = (X(P_1)Z(P_1) - Z(P_1)X(P_2)), $$

by the differential addition formulas from above. Lastly,

$$ \kappa_{01} = 2((X(P_1)X(P_2) + Z(P_1)Z(P_2))(X(P_1)Z(P_2) + X(P_2)Z(P_1)) + 2AX(P_1)X(P_2)Z(P_1)Z(P_2)) $$

by homogenisation of [BDLS20, Example 4.4].

For the three way addition, Section 4.9.4 already gives the formula for Z. We can obtain $X(P_1 + \tilde{P}_2 + Q)$ from the formula for Z simply by replacing Q by $Q + T$:

$$ X(P_1 + \tilde{P}_2 + Q) = \frac{(X(\tilde{P}_2)X(P_1 + Q) - Z(P_1 + Q)Z(\tilde{P}_2))(X(\tilde{P}_1)X(P_2 + Q) - Z(P_2 + Q)Z(\tilde{P}_1))}{(X(P_1 - P_2)X(\tilde{Q}) - Z(P_1 - P_2)Z(\tilde{Q}))} $$

Algorithms 5.7 and 5.8 give the algorithms for the compatible addition and the three way cubical addition in the Montgomery model respectively.

We can combine Algorithms 5.7 and 5.8 to compute the addition step of Algorithm 4.7. We compute $P_1 + P_2$ by a compatible addition, then take an arbitrary lift $P_1 + P_2$ and then apply the three way addition to obtain $P_1 + \tilde{P}_2 + Q$. We remark that in Algorithm 5.8, we already have $P_1 - \tilde{P}_2$ from Algorithm 5.4. Furthermore, while we describe these algorithms with divisions for simplicity, we can use Lemma 4.15 to remove all divisions: we can multiply $P_1 + \tilde{P}_2, P_1 + \tilde{P}_2 + Q$ by the same projective factor $\lambda$ and still get the same biextension element $g_{P_1+P_2,Q}$.
**Input:** $\tilde{P}_1 = (X(P_1), Z(P_1)), \tilde{P}_2 = (X(P_2), Z(P_2)), P_1 \tilde{+} Q = (X(P_1+Q), Z(P_1+Q)), P_2 \tilde{+} Q = (X(P_2+Q) : Z(P_2+Q)), P_1 \tilde{+} P_2 = (X(P_1+P_2), Z(P_1+P_2))$

**Output:** $P_1 + \tilde{P}_2 + Q = (X(P_1 + P_2 + Q), Z(P_1 + P_2 + Q)) = \text{ThreeWayAdd}(\tilde{P}_1, \tilde{P}_2, \tilde{Q}, P_1 \tilde{+} Q, P_2 \tilde{+} Q, P_1 \tilde{+} P_2)$

→ Compute $(X(P_1 \tilde{-} P_2), Z(P_1 \tilde{-} P_2)) = \text{DiffAdd}(\tilde{P}_1, \tilde{P}_2, P_1 \tilde{+} P_2)$

→ $X(P_1 + \tilde{P}_2 + Q) = \frac{(X(P_2)X(P_1+Q) - Z(P_1+Q)Z(P_2))(X(P_1)X(P_2+Q) - Z(P_2+Q)Z(P_1))}{X(P_1-P_2)X(Q) - Z(P_1-P_2)Z(Q)}$

→ $Z(P_1 + \tilde{P}_2 + Q) = \frac{(X(P_2)Z(P_1+Q) - X(P_1+Q)Z(P_2))(X(P_1)Z(P_2+Q) - X(P_2+Q)Z(P_1))}{X(P_1-P_2)X(Q) - Z(P_1-P_2)Z(Q)}$

<caption>ALGORITHM 5.8. Three way cubical additions in the Montgomery model</caption>

In practice, this means we compute the numerators $n$ and denominators $d$ of $X(P_1 + \tilde{P}_2 + Q)$ and $Z(P_1 + \tilde{P}_2 + Q)$ separately, and use the formulas

$$X(P_1 + P_2) = \kappa_{00}Z(P_1 - P_2)d(X(P_1 + P_2 + Q))d(Z(P_1 + P_2 + Q)),$$

$$Z(P_1 + P_2) = \kappa_{11}X(P_1 - P_2)d(X(P_1 + P_2 + Q))d(Z(P_1 + P_2 + Q)),$$

$$X(P_1 + P_2 + Q) = n(X(P_1 + P_2 + Q))d(Z(P_1 + P_2 + Q))X(P_1 - P_2)Z(P_1 - P_2),$$

$$Z(P_1 + P_2 + Q) = n(Z(P_1 + P_2 + Q))d(X(P_1 + P_2 + Q))X(P_1 - P_2)Z(P_1 - P_2).$$

We could use a similar strategy to compute a biextension double $2\tilde{P}, 2\tilde{P} + Q$, but it is faster to use a cubical double and differential addition for this case, as explained in Section 4.7.

These algorithms combined give us the $32M + 4S + 2m_0$ count for an addition. There is probably room for improvement for these formulas by rearranging the arithmetic operations. Also it could potentially be faster to do a compatible addition to compute $P_1 + P_2 + Q$ projectively rather than $P_1 - P_2$, and then using cubical arithmetic to get $P_1 + \tilde{P}_2$ from a choice of $P_1 + \tilde{P}_2 + Q$.

We refer to the implementation in [Rob23b] for a variant which combines a double and addition into a DoubleAndAdd, which costs $17M + 8S + 3m_0$ (so strangely, is faster than just an addition, but this is probably an artefact of the fact that the addition is not at all optimised).

5.2.3. Results. Giacomo Pope kindly implemented our algorithm in Rust. Section 5.2.3 give some timings on a Intel Core i5-1335U, with turboboost disabled, on a supersingular curve $E$ over $\mathbb{F}_{p^2}$ with $p = 2^{604} \cdot 3^{363} - 1$ and $p = 2^{74} \cdot 3^{41} - 1$ respectively.

<table><thead><tr><th></th><th>2<sup>604</sup></th><th>3<sup>363</sup></th><th>2<sup>604</sup> · 3<sup>363</sup></th><th>2<sup>74</sup></th><th>3<sup>41</sup></th><th>2<sup>74</sup> · 3<sup>41</sup></th></tr></thead><tbody><tr><td>Tate pairing</td><td>26.6 ms</td><td>34.7 ms</td><td>57.9 ms</td><td>80.6 µs</td><td>100.8 µs</td><td>166.4 µs</td></tr><tr><td>Weil pairing</td><td>35.2 ms</td><td>50.9 ms</td><td>103.84 ms</td><td>112.8 µs</td><td>149.7 µs</td><td>300.9 µs</td></tr></tbody></table>

One nice advantage of Theorem 1.1 is that the cubical ladder, since it is derived from the Montgomery ladder, is naturally constant time. Constant time constraints are different for pairings than for standard elliptic curve cryptography: in elliptic curve cryptography the scalar is secret, while the base point may or may not be public. For pairings, the scalar is usually public, but the paired points may be secret. The pairing arithmetic such need to not
depend on the paired point. By Remark 5.3, this is naturally the case for the cubical ladder in the Montgomery model, except for a few exceptions (which do not happen when $\ell$ is odd).

Compared to a previous constant time implementation using Miller's algorithm, our new algorithm is between 2.5× to 3× faster over $\mathbb{F}_{p^2}$ for $p = 2^{74} \cdot 3^{41} - 1$.

A last advantage of Theorem 1.1 is that in isogeny based cryptography, we often represent a basis $(P_1, P_2)$ of $N$-torsion by $x(P_1), x(P_2), x(P_1 + P_2)$ (or via their $(X:Z)$ projective coordinates). Indeed, it usually does not matter if we compute an isogeny $\phi$ or its opposite $-\phi$, so it is harmless to replace the basis $(P_1, P_2)$ by $(-P_1, -P_2)$, hence represent it by $x(P_1), x(P_2), x(P_1 + P_2)$. Since Theorem 1.1 naturally take such a data as input, it allows to compute the pairing $e(P_1, P_2)$ without the need of taking a square root first.

**5.3. The theta model.** Explicit formulas in the level 2 theta model (in any dimension) are given in [LR16]. Although it was not stated in this form there, the formulas compute the correct cubical arithmetic.

In the theta model, a cubical differential addition (with the base point normalised) costs $2S + 3M$, and a doubling costs $3S + 1M + 2m_0$. Since the cubical biextension ladder uses one doubling and two differential additions by bits, the total cost is $7S + 7M + 2m_0$.

There is a variant which costs $2S + 3M + 1m_0$ for a differential addition, and $4S + 2m_0$ for a doubling; using this variant the cubical biextension ladder than costs $8S + 6M + 4m_0$ by bits.

**5.4. Short Weierstrass curves.** On a short Weierstrass curve $y^2 = x^3 + ax + b$, $x$-only doublings and differential additions formulas are given by:

$$x(P+Q)x(P-Q) = \frac{-4b(x(P) + x(Q)) + (x(P)x(Q) - a)^2}{(x(P) - x(Q))^2} = \frac{-4bZ(P)Z(Q)(X(P)Z(Q) + X(Q)Z(P)) + (X(P)X(Q))^2}{(X(P)Z(Q) - X(Q)Z(P))^2}$$

$$x(2P) = \frac{(x(P)^2 - a)^2 - 8bx(P)}{4(x(P)^3 + ax(P) + b)} = \frac{(X(P)^2 - aZ(P))^2 - 8bX(P)Z(P)^3}{4Z(P)(X(P)^3 + aX(P)Z(P)^2 + bZ(P)^3)}$$

Exactly as for Montgomery curves the denominators are precisely the ones coming from the cubical arithmetic:

$$Z(P+Q)Z(P-Q) = (X(Q)Z(P) - X(P)Z(Q))^2, \\ Z(2P) = 4(X(P))^3Z(P) + aX(P)Z(P)^3 + bZ(P)^4.$$

This means that the standard ladder [B]02 on short Weierstrass curves already compute the correct cubical arithmetic (modulo the fact that, in the ladder, we need to divide the $X$ coordinate by $x(P_0)$ where $P_0 = (x(P_0), 1)$ is the base point, rather than multiplying the $Z$ coordinate by $x(P_0)$, in order to keep track of the correct cubical factor).

Thus a cubical doubling costs $3M + 4S + 2m_0$ and a cubical differential addition costs $6M + 2S + 2m_0$, so a cubical biextension ladder costs a quite expensive $15M + 8S + 6m_0$ by bits. This allows to compute pairings on short Weierstrass curves knowing only $x(P), x(Q), x(P+Q)$.

We leave for future work improving these formulas. Maybe keeping track of the $Y$ coordinate could help. An alternative approach would be to use a Montgomery ladder in co-Z coordinates [Mel07; GJMRV11; Ham20], since these gives faster formula. But we would need to extend this co-Z approach to a 3 point ladder (i.e., a ladder computing $\ell P, \ell P + Q$), and we would still need to track the correct cubical factor for each point somehow (scaling to the same co-Z coordinates on $mP, (m+1)P, (m+1)P + Q$ kills the cubical information if we do not keep it somewhere else).
## 6. APPLICATIONS

We discuss various applications of the cubical arithmetic. In Section 6.1 we look at pairing based cryptography and how the cubical arithmetic on Montgomery models from Section 5 compares with the usual Miller's algorithm. We also derive new (to my knowledge), more or less interesting, formulas for the standard Miller algorithm by using the cubical point of view.

In Section 6.2, we briefly mention some applications of cubical arithmetic beyond pairings, notably on isogenies and radical isogenies.

In Section 6.3, we reinterpret, using cubical arithmetic, Doliskani's supersingularity test (see [Dol18; BGS22]) as a self Tate pairing test.

Finally, our most important application is probably Section 6.4 where we use cubical arithmetic to obtain a new powerful side channel attack against the Montgomery ladder for Montgomery curves. Namely, one projective coordinate leak in the Montgomery ladder allows to solve the elliptic curve dlp by reduction to the base field dlp.

Each of these applications would deserve an article on its own<fnref n="2" />. We won't detail these applications much in this paper, because it is quite long already.

### 6.1. Pairing based cryptography.

6.1.1. *Comparison with Miller's algorithm for pairing based cryptography.* For pairing based cryptography on elliptic curves, with embedding degree $d > 1$, it is convenient to use the Tate pairing with $P \in \mathbb{G}_1 \subset E(\mathbb{F}_q)$, $Q \in \mathbb{G}_2 \subset E(\mathbb{F}_{q^d})$, and $d$ even to allow for denominator elimination. We recall that $Q \in \mathbb{G}_2$ if and only if $\pi_q(Q) = qQ$.

Counting only operations involving the big field $\mathbb{F}_{q^d}$, Miller's algorithm, with denominator elimination, costs $1M + 1S + 1m$ by doubling, and $1M + 1m$ by addition. Here $1m$ denotes a "mixed multiplication", meaning a multiplication between a coefficient in $\mathbb{F}_q$ and a coefficient in $\mathbb{F}_{q^d}$.

When denominator elimination is not possible (because $d$ is odd or $Q$ is not in $\mathbb{G}_2$), the cost becomes $2M + 2S + 1m$ by doubling, and $2M + 1m$ by addition. (In practice, working with projective coordinates for $P$, the number of mixed multiplications $m$ is a bit higher, see Section 6.1.2.) On the other hand, when $d$ is odd, there are variants of Miller's algorithm in the literature (notably [BELL10]) which we will discuss in Section 6.1.2, which achieve $2S + 1M + 3m$ for a doubling and $1M + 2.5m$ for an addition.

Using our arithmetic of biextension on Kummer lines in the Montgomery model, only counting the operations on the big field, we have $2S + 1M + 2m$ by bit for the cubical biextension ladder. This is better than Miller's algorithm (even the improved variant), except when denominator elimination is available.

Our main difficulty is that it is not clear how to do an efficient denominator elimination in cubical arithmetic when the embedding degree is even. If $d$ is even and $Q \in \mathbb{G}_2$, then $x(Q)$ lies in a strict subfield of $\mathbb{F}_{q^d}$. This allows to do denominator elimination, i.e. not compute the denominator of the Miller functions $\mu_{uP,vP}$ evaluated at $Q$, since this evaluated denominator is $x(Q) - x((u + v)P)$ is in a strict subfield, hence killed by the final exponentiation.

In Theorem 1.1, for $P \in \mathbb{G}_1$, $Q \in \mathbb{G}_2$ and $d$ even, then while $x(P), x(Q)$ lies in a subfield, we also need to use $x(P + Q)$ which does not. This means that for the cubical biextension

<footnote n="2">²Which I cannot promise I will ever write someday; I already made such promises in previous papers and so far I have a bad track record on keeping these... The only exception will hopefully be the monodromy leak attack, because several people have encouraged me to publish it in its own paper, rather than to hide it at the end of a paper about biextensions.</footnote>
ladder, the cost is $2S + 1M + 2m$ when the current bit is 1 and the difference point is $P + Q$, and $2S + 1M' + 2m$ when the current bit is 0 and the difference point is $Q$, where $M'$ means that we multiply by a coordinate in $\mathbb{F}_{p^{d/2}}$. So even embedding degree $d$ still helps for the cubical ladder.

Another way it helps, which is well known, is that an inverse in $\mathbb{F}_{q^d}$ can be written $1/x = \bar{x}/(x\bar{x})$ where $\bar{x}$ is the quadratic Galois conjugate of $x$ and in particular $x\bar{x} \in \mathbb{F}_{q^d/2}$. In situations where $x\bar{x}$ will be killed by the final exponentiation, this allows to use $\bar{x}$ instead of $1/x$, hence replace a division by (one or several) multiplications.

It could be interesting to work out formulas where we represent $m\widetilde{P} + Q$ using only $Z(m\widetilde{P} + Q)$ rather than via $(X(m\widetilde{P} + Q), Z(m\widetilde{P} + Q))$.

On the other hand, for the Ate or optimal Ate pairing, since the computation is done with $P \in \mathbb{G}_2$ and $Q \in \mathbb{G}_1$, it is plausible that the cubical arithmetic could be faster than the usual Miller's algorithm, since the overall number of operations for a generic pairing is much smaller.

### 6.1.2. Cubical arithmetic and new formulas for Miller's algorithm.
In this section, rather than comparing Miller's algorithm with the cubical arithmetic on Kummer lines, we will compare Miller's algorithm with the cubical arithmetic on elliptic curves, where we represent $\widetilde{P}$ by using the level 3 projective Weierstrass coordinates $P = (X(P) : Y(P) : Z(P))$ along with the level 1 affine cubical coordinate $Z_1(\widetilde{P})$.

To simplify, we will consider the case of the Tate pairing $P \in \mathbb{G}_1, Q \in \mathbb{G}_2$. In Miller's algorithm, we compute $f_{\ell,P}((Q) - (0))$. In the cubical arithmetic, we compute $\ell\widetilde{P}, \ell\widetilde{P} + Q$. The relationship between the two approaches is given by Porism 4.25. In particular, assuming we fix $\widetilde{P}, \widetilde{Q}, \widetilde{P} + Q$ so that $Z_1(\widetilde{P}) = Z_1(\widetilde{Q}) = Z_1(\widetilde{P} + Q) = 1$, under our assumptions on $P, Q$, we have that $Z_1(\ell\widetilde{P} + Q)$ is equal to $f_{\ell,P}((Q) - (0))$ up to a factor lying in the small field $\mathbb{F}_q$, which will be killed by the final exponentiation.

Now we already saw in Remark 3.16 that the biextension arithmetic provided new insight on how to compute $f_{\ell,P}$ via the functions $\mu_{uP,vP}$. We will now derive new formulas using the cubical arithmetic.

Let's start by formulas to compute $Z_1(\ell\widetilde{P} + Q)$, using Section 4.9.3. We have the doubling formula

$$Z_1(2m\widetilde{P}) = 2y(mP)^2 Z_1^4(m\widetilde{P})$$

and the differential addition formulas

$$ \begin{aligned} Z_1((2m\widetilde{P}+1)P)Z_1(\widetilde{P}) &= Z_1((m\widetilde{P}+1)P)^2 Z_1(m\widetilde{P})^2 (x(mP) - x((m+1)P)), \\ Z_1(2m\widetilde{P}+Q)Z_1(\widetilde{Q}) &= Z_1(m\widetilde{P}+Q)^2 Z_1(m\widetilde{P})^2 (x(mP) - x(mP+Q)), \\ Z_1((2m\widetilde{P}+1)P+Q)Z_1(\widetilde{P}+Q) &= Z_1((m\widetilde{P}+1)P+Q)^2 Z_1(m\widetilde{P})^2 (x(mP) - x((m+1)P+Q)). \end{aligned} $$

This allows to compute $\ell\widetilde{P} + Q$ via a cubical ladder using the points $m\widetilde{P}, (m+1)P, (m+1)\widetilde{P} + Q$. (Note the similarity with Remark 3.16).

We could also use a double and add algorithm, using

$$Z_1(2m\widetilde{P}+Q)Z_1(\widetilde{Q}) = Z_1(m\widetilde{P}+Q)^2 Z_1(m\widetilde{P})^2 (x(mP) - x(mP+Q))$$

for doublings, and the three way addition

$$Z_1(P + m\widetilde{P} + Q) = \frac{Z_1(P + m\widetilde{P})Z_1(m\widetilde{P} + Q)Z_1(P + Q)}{Z_1(\widetilde{P})Z_1(m\widetilde{P})Z_1(\widetilde{Q})} \frac{l_{mP,P}(Q)}{(x(Q) - x(mP))(x(Q) - x(P))'}$$
for an addition. We could also compute the doubling using the following three way addition:

$$Z_1(mP + \widetilde{mP} + Q) = \frac{Z_1(m\widetilde{P} + mP)Z_1(m\widetilde{P} + Q)^2}{Z_1(\widetilde{mP})^2Z_1(\widetilde{Q})} \frac{l_{mP,mP}(Q)}{(x(Q) - x(mP))^2}$$

I implemented both methods and the pairing is computed correctly, but these formulas seem a priori no faster than Miller's standard algorithm.

However, the following approach seems the most promising: use a double and add algorithm, using

$$Z_1(2m\widetilde{P} + Q)Z_1(\widetilde{Q}) = Z_1(m\widetilde{P} + Q)^2Z_1(\widetilde{mP})^2(x(mP) - x(mP + Q))$$

as before for doublings, and

$$Z_1((2m + \overline{1})P + Q)Z_1(P + \widetilde{Q}) = Z_1((m + \overline{1})P + Q)^2Z_1(\widetilde{mP})^2(x(mP) - x((m+1)P + Q))$$

for a doubling and addition.

Now we remark that the function $Q \mapsto x(mP) - x(mP + Q)$ has for divisor $(0_E) + (-2mP) - 2(-mP)$ so is equal to a biextension function $g_{mP,mP}(Q)$. In particular, $(x(mP) - x((m+1)P + Q))$ is given by $g_{mP,mP}(P + Q)$.

This gives the following alternative strategy to compute the (normalised) Miller function $f_{\ell,P}(Q)$. Usually we use $f_{u+v,P}(Q) = f_{u,P}(Q)f_{v,P}(Q)\mu_{uP,vP}(Q)$ with $\mu_{uP,vP}$ the normalised function with divisor $D_{uP+vP} - D_{uP} + D_{vP}$. In particular, the above formula gives the same Miller doubling: $f_{2u,P}(Q) = f_{u,P}(Q)^2\mu_{uP,uP}(Q)$. But we also have the following Miller DoubleAndAdd: $f_{2u+1,P}(Q) = C f_{u+1,P}(Q)^2\mu_{uP,uP}(P + Q)$. Here the constant C refers to the fact that the function $Q \mapsto \mu_{uP,uP}(P + Q)$ is no longer normalised; however under our assumptions that $P \in \mathbb{G}_1$, C is in a strict subfield so will be killed by the final exponentiation. Of course, this formula can be proved directly without going through the cubical arithmetic: the left hand side has divisor $D_{(2u+1)P} - (2u+1)D_P$ while the right hand side has divisor $2D_{(u+1)P} - 2(u+1)D_P + t_P^*(D_{2uP} - 2D_{uP}) = 2D_{(u+1)P} - 2(u+1)D_P + D_{(2u+1)P} - 2D_{(u+1)P} + D_P = D_{(2u+1)P} - (2u+1)D_P$.

Using the standard sign conventions, the latter formulas reads:

$$f_{2u+1,P}(Q) = C f_{u+1,P}(Q)^2 \mu_{uP,uP}(Q - P).$$

This gives the following double and DoubleAndAdd algorithm: at each step we compute $mP$ or $(m+1)P$ along with $f_{m+1,P}(Q)$ up to some factor in a subfield.

We use either: $f_{2m+1,P}(Q) = C_1 f_{m+1,P}(Q)^2 \mu_{mP,mP}(Q - P)$ (in which case we need to compute $mP$ if we had $(m+1)P$, and we also compute $2mP$) or $f_{2m+2,P}(Q) = C_2 f_{m+1,P}(Q)^2 \mu_{(m+1)P,(m+1)P}(Q)$ (in which case we need to compute $(m+1)P$ if we had $mP$, and we also compute $2(m+1)P$). This makes a DoubleAndAdd very similar to a simple doubling.

We compute in this way $f_{\ell,P}(Q)$ (or $f_{\ell+1,P}(Q)$ if we prefer). As far as I know, this DoubleAndAdd method is new. We refer to [Rob23b] for the implementation.

Let us compare this strategy with known optimisation of Miller's algorithm. In [BELL10], they use the formula $f_{u+v,P}(Q) = \frac{1}{f_{-u,P}(Q)f_{-v,P}(Q)l_{-uP,-vP}(Q)}$ Compared to the usual formula which involves $\mu_{uP,vP}(Q) = l_{uP,vP}(Q)/v_{uP+vP}(Q)$, where $v_{uP+vP} = l_{uP+vP,-uP-vP}$, this saves the denominator $v_{uP+vP}(Q)$ of $\mu_{uP,vP}$.

In [DZZZ22] the authors introduce functions with divisors $\ell(P) + (-\ell P) - (\ell + 1)(0)$ which give a streamlined double and add formula. We remark that such a function is given by $f_{-\ell,P}$, so [DZZZ22] is subsumed by the formula above from [BELL10].
Under the assumption $P \in \mathbb{G}_1, Q \in \mathbb{G}_2$ and counting only operations in the big field, [BELL10] have a complexity (for $d$ odd) of $2S + 1M + 3m$ for a doubling (compared to $2S + 2M + 5m$ for the classical Miller formula), and $1M + 2.5m$ for an addition (compared to $2M + 5m$ for the classical Miller formula); so a DoubleAndAdd costs $2M + 2S + 5.5m$ (compared to $4M + 2S + 10m$). When $d$ is even, they have the same complexity (in the big field) as the usual Miller's formula of $1M + 1S + 1m$ for a doubling and $1M + 1m$ for an addition; so a DoubleAndAdd costs $2M + 1S + 2m$. (We assume here that multiplying an element in $\mathbb{F}_p$ with an element of $\mathbb{F}_{p^d/2}$ is half the cost of a multiplication of an element in $\mathbb{F}_p$ with an element in $\mathbb{F}_{p^d}$.)

Our new DoubleAndAdd method cost $2S + 2M + 5m$ both for a doubling and a DoubleAndAdd. When $d$ is even, we have the same difficulty as in Section 6.1.1 that $x(P+Q)$ does not lies in a subfield. If $\mathbb{F}_{p^d} = \mathbb{F}_{p^{d/2}}[i]$ with $i^2 = c$, we can use that $1/(x + iy) = (x - iy) / (x^2 - dy^2)$ where $x, y \in \mathbb{F}_{p^d/2}$ and since $x^2 - dy^2$ is in a strict subfield, this allows to replace divisions by multiplications. The cost becomes $2M + 1S + 3m$ for a DoubleAndAdd which involves $x(P + Q)$ and $1M + 1S + 1m$ for a doubling which involves $x(Q)$ which is in a subfield. So for $d$ odd, the DoubleAndAdd method is better than the usual Miller's formula but not better than [BELL10], and for $d$ even roughly similar to [BELL10].

However we can combine the DoubleAndAdd variant with the following line folding trick (see [LL11, Lemma 3.2] and the references): $\frac{l_{T,T}(Q)}{v_l^2(Q)v_{2T}(Q)} = \frac{1}{l_{-T,-T}(Q)}$. This allows to delay a vertical line evaluation and fold it in a standard line evaluation. With the DoubleAndAdd algorithm, since we change the evaluated point according to whether the current bit is 0 or 1, we can only apply this line folding trick for consecutive 0s or 1s. Using this line folding trick, [LL11] obtain a complexity, for $d$ odd and forgetting mixed additions, of $2S + 1M$ for a doubling and $2S + 2M$ for a double and add. (They give $2S + 1M$ or $2S + 2M$ according to which branch of the algorithm is evaluated, but it seems to me that in the branch where a parabola is computed one should account that evaluating this parabola costs $2M$). By comparison, our DoubleAndAdd algorithm (and also not counting mixed additions), we have a cost of $2S + 1M$ for consecutive doublings or DoubleAndAdds, and, as we have seen above, a cost of $2S + 2M$ when we switch from a bit equal to 0 to a bit equal to 1 or conversely. In the case when $d$ is even, using the line folding trick, consecutive DoubleAndAdds cost $1M+1S$, the same as doublings (whether they are consecutive or not). So this DoubleAndAdd algorithm could be interesting when $l$ has low Hamming weight.

**6.2. Algorithms for abelian varieties.** The cubical arithmetic we develop has important applications for algorithms on abelian varieties beside pairings.

Notably, we will see that explicit formulas for the cubical arithmetic gives explicit formulas for the theta group action (at level $l$). Using the generic isogeny framework from [Rob21a, § 2.9, § 4.2], this explicit action gives formulas for isogenies (or even allows to compute a basis of algebraic theta functions or to change levels). We have also just seen in this paper how the cubical arithmetic allows to compute pairing.

Now by the construction of the algebraic Riemann relations from Proposition 4.1, the cubical arithmetic can be derived from explicit formulas for the theorem of the square: $t_{P+Q}^*\mathcal{L} \otimes \mathcal{L} \simeq t_P^*\mathcal{L} \otimes t_Q^*\mathcal{L}$, which also (implicitly) encodes addition formulas. In summary: explicit formulas for the theorem of the square, which induces explicit formulas for the cubical arithmetic, is key for all standard algorithms on abelian varieties: additions, pairings, isogenies.
6.2.1. **Theta group arithmetic of higher level through cubical arithmetic.** First, we show that cubical arithmetic gives a convenient way to work out the arithmetic of theta group of level $ℓn$ while using coordinates of level $n$ evaluated on cubical points of $ℓ$-torsion. As explained in [LR22a], having this explicit theta group action can then be used to compute isogenies by taking suitable traces; so the cubical arithmetic can be used to find isogeny formulas in any models (by contrast to [LR12; CR15; LR22a] which use the theta model).

The basic idea is as follows: let $(A, \mathcal{L})$ be a polarised abelian variety. To compute $ℓ$-isogenies, we need to work out the arithmetic of the theta group $G(\mathcal{L}^\ell)$. We recall from Section 2 that an element $(P, g_{\mathcal{L}^\ell, P}) \in G(\mathcal{L}^\ell)$ is such that $P \in \text{Ker } \Phi_{\mathcal{L}^\ell}$ where $\Phi_{\mathcal{L}^\ell}(P) = t_P^* \mathcal{L}^\ell \otimes \mathcal{L}^{-\ell}$ is the polarisation associated to $\mathcal{L}^\ell$ and $g_{\mathcal{L}^\ell, P}$ is a section of $t_P^* \mathcal{L}^\ell \otimes \mathcal{L}^{-\ell}$.

In particular, a cubical point $\tilde{P}$ for $\mathcal{L}^\ell$ is enough to determine $g_{\mathcal{L}^\ell, P}$. Now a cubical point $\tilde{P}$ for $\mathcal{L}$ does determine a cubical point for $\mathcal{L}^\ell$, simply by taking the tensor product to the $\ell$ of the rigidification at $P$; we will denote this by $\tilde{P}^{\otimes \ell}$. (In particular, $\zeta \tilde{P}$ induce the same cubical point for $\mathcal{L}^\ell$ whenever $\zeta^\ell = 1$.) We remark also that if $[\tilde{P}_1, \tilde{P}_2, \tilde{P}_3, \tilde{P}_4; \tilde{Q}_1, \tilde{Q}_2, \tilde{Q}_3, \tilde{Q}_4]$ are in Riemann position, then certainly their $\ell$-fold tensor product are in Riemann position; so $\tilde{P} \mapsto \tilde{P}^{\otimes \ell}$ is compatible with the cubical arithmetic with respect to $\mathcal{L}$ and $\mathcal{L}^\ell$ respectively. Furthermore, if $\gamma$ is the function associated to the Riemann relation $[P_1, P_2, P_3, P_4; Q_1, Q_2, Q_3, Q_4]$ for $\mathcal{L}$, then $\gamma^\ell$ is the corresponding function for $\mathcal{L}^\ell$. Finally, if $X_1, \dots, X_\ell$ are sections of $\Gamma(\mathcal{L})$, then $X = \prod_{i=1}^\ell X_i$ is a section of $\Gamma(\mathcal{L}^\ell)$, and $X(\tilde{P}^{\otimes \ell}) = \prod_{i=1}^\ell X_i(\tilde{P})$.

This allows to work on the theta group $G(\mathcal{L}^\ell)$ (i.e. at level $\ell$) while using cubical points for $\mathcal{L}$ (i.e. at the base level). (We remark that the method of Remark 2.11, to compute the $\ell m$ Tate pairing for $D_1$ when working with $D = mD_1$ goes in the other direction: we work with cubical points for $D$, but we try to normalise things such that the cubical points come from the $m$-fold tensor product of cubical points on $D_1$; that's how I found out this trick!)

We illustrate this approach in the following situation. Let $K \subset A[\ell]$ be an isotropic subgroup. By [Rob21a, Chapter 2], to compute isogeny formulas, we first need to compute a lift $\tilde{K}$ of $K$ to the theta group $G(\mathcal{L}^\ell)$, and then compute the action of $\tilde{K}$ on sections of $\Gamma(\mathcal{L}^\ell)$.

If $P \in A[\ell]$, and $\tilde{P}$ is a cubical point with respect to $\mathcal{L}$ then we have seen above that $\tilde{P}^{\otimes \ell}$ is a cubical point with respect to $\mathcal{L}^\ell$, hence encodes a theta group element $g_P \in G(\mathcal{L}^\ell)$. This element $g_P$ acts on the sections of $\Gamma(\mathcal{L}^\ell)$. Here is how we can recover this action using the cubical point $\tilde{P}$, in the particular case when the section $u \in \Gamma(\mathcal{L}^\ell)$ can be written as $u = \prod_{i=1}^\ell u_i$, where $u_i \in \Gamma(\mathcal{L})$.

First, we need to take $\tilde{P}$ such that $g_P$ is an element of $\ell$-torsion. If $\mathcal{L}$ is symmetric, it if often convenient to add the further constraint that $g_P$ should be symmetric: when $\mathcal{L}$ is symmetric there is a canonical involution $\delta_{-1}$ on $G(\mathcal{L}^\ell)$ defined in [Mum66], and an element is symmetric if $\delta_{-1}g_P = g_P^{-1}$. We will see below the definition of a canonical lift of $\ell$-torsion $\tilde{P}$, and that such a $\tilde{P}$ gives an element $\tilde{P}^{\otimes \ell}$ such that the associated theta group element $g_P \in G(\mathcal{L}^\ell)$ is symmetric of $\ell$-torsion.

Next, for a $R \in A$, we take an arbitrary cubical point $\tilde{R}$ for $\mathcal{L}$, and we take $R + P$ such that $R + \ell P = \tilde{R}$. By homogeneity, this equation constrains $R + P$ from $\tilde{R}, \tilde{P}$ completely, up to the action by $\mu_\ell \subset \mathbb{G}_m$.

Now the theta group action of $g_P$ on $u = \prod u_i$ is induced by the value $\prod u_i(\tilde{P} + \tilde{R})$; since we take an $\ell$-fold product the action of $\mu_\ell$ above is killed, hence this product is well defined. More precisely, if $X_0$ has for divisor $D$, then $g_P \in G(\ell D)$ can be described by the formula $g_P(R) = X_0(\tilde{R} + P)^m / X_0(\tilde{R})^m$.
6.2.2. Canonical lifts of $\ell$-torsion. As a warm up, we first look at the theta group arithmetic on $G(D)$ encoded by cubical points $\tilde{P}$ for $P \in A[D]$. Then we will extend this to the theta group arithmetic for $G(\ell D)$ encoded by $\tilde{P}^{\otimes \ell}$ for $P \in A[\ell]$.

In this section, we freely use the notations $\delta_{n}, e_{*}, \dots$ from [Mum66]. Let's assume that we are on level $n$ for simplicity, so $P$ is a point of $n$-torsion, and we assume further that $D$ is symmetric, and if $n$ is even that (the line bundle associated to) $D$ is totally symmetric. In particular, $e_{D, *}(T) = 1$ for all $T \in A[2]$, so the symmetric elements $g_T \in G(D)$ above $T$ are of order 2. Furthermore, over $\bar{k}$, the theta group $G(D)$ is isomorphic to the Heisenberg group $\bar{k}^* \times (\mathbb{Z}/n\mathbb{Z})^g \times \mu_n^g$. On the Heisenberg group, we have that $\delta_n(\alpha, x, \zeta) = (\alpha^{n^2}, nx, \zeta^n)$, and that $(\pm 1, x, \zeta)$ are the two symmetric elements above $(x, \zeta)$.

One need to be careful that there are two different laws for the cubical arithmetic when $P \in A[D]$: the first one comes the action by $G(D)$. Iterating this action $m$ times we get the cubical point $\tilde{P} \cdot \tilde{P} \cdot \dots \tilde{P}$. If $\tilde{P}$ corresponds to $g_P \in G(D)$, this laws corresponds, by definition of the cubical action, to the theta group law $g_P^m$. The second law comes from the cubical exponentiation: $m\tilde{P}$. One can check (the easiest way is to use the analytic interpretation of the cubical law through the canonical factor of automorphy, see [LR22a]) that $m\tilde{P}$ corresponds to $\delta_m g_P$. (A way to reformulate this is as follows: in the cubical ladder, if $\tilde{P}$ corresponds to $g_P \cdot \tilde{\emptyset}$, then $-\tilde{P}$ corresponds to $\delta_{-1}g_P \cdot \tilde{\emptyset}$. That's why the ladder will give $n\tilde{P} = \delta_n g_P \cdot \tilde{\emptyset}$. But the cubical arithmetic is also compatible with the group action by Lemma 4.12. So in the ladder if we would use $\tilde{P}^{-1} := g_P^{-1} \cdot \tilde{\emptyset}$ instead of $-\tilde{P}$, we would obtain $\tilde{P}^n = g_P^n \cdot \tilde{\emptyset}$.)

In particular, $n\tilde{P} = \tilde{\emptyset}$ corresponds to $\delta_n g_P = 1$ which is equivalent to $g_P$ being of $n^2$-torsion. This is why the equation $n\tilde{P} = \tilde{\emptyset}$ has $n^2$-solutions over $\bar{k}$.

If $n = 2n' + 1$ is odd, we can check that $g_P$ is of $n$-torsion is equivalent to $\delta_{n'+1}g_P = \delta_{-n'}g_P$. This corresponds to the equation $(n' + 1)\tilde{P} = -n'\tilde{P}$, which indeed has $n$ solutions over $\bar{k}$. If $n = 2n' + 2$ is even, we have that $g_P$ is of $2n$-torsion is equivalent to $\delta_{n'+2}g_P = \delta_{-n'}g_P$, which corresponds to the equation $(n' + 2)\tilde{P} = n'\tilde{P}$ and has $2n$ solutions over $\bar{k}$.

Now let's consider a point $P \in A[\ell]$, and take a cubical point $\tilde{P}$. By Section 4.2.5, we have if $\ell \tilde{P} = \lambda_0 \tilde{\emptyset}$, $(\ell + 1) \tilde{P} = \lambda_0 \lambda_1' \tilde{P}$, with $\lambda_0^2 = \lambda_1'^{\ell}$.

Let us first assume that $\ell = 2\ell' + 1$ is odd. Then there is a unique $\beta$ such that $\beta^\ell = \lambda_0$, $\beta^2 = \lambda_1'$, $\beta$ is rational since $\lambda_0, \lambda_1'$ are rational, and we have $(u\ell + v)\tilde{P} = \beta^{\ell u^2+2v}\tilde{P}$. Furthermore, it is easy to check that replacing $\tilde{P}$ by $\gamma\tilde{P}$ changes $\beta$ by $\beta\gamma^\ell$, hence $\lambda_0$ by $\lambda_0\gamma^{\ell^2}$ and $\lambda_1'$ by $\lambda_1'\gamma^{2\ell}$. Also, by the monodromy interpretation of the Tate pairing, $\lambda_1' = \beta^2$ is (a representative of the class of) the non reduced Tate pairing $e_{T,\ell}(P, P)$.

From this, we see that the equation $\ell \tilde{P} = \tilde{\emptyset}$ (which has $\ell^2$ solutions over $\bar{k}$) imposes $\beta^\ell = 1$, which is not enough to guarantee that $(\ell + 1)\tilde{P} = \tilde{P}$. However, both equations $\ell \tilde{P} = \tilde{\emptyset}$, $(\ell + 1)\tilde{P} = \tilde{P}$ are enough to impose $\lambda_0 = \lambda_1' = 1$, which imposes $\beta = 1$ in the odd case, and $\beta^2 = 1$ in the even case. If these equations are satisfied, we say that $\tilde{P}$ is a canonical cubical point of $\ell$-torsion; we will also say it is a canonical lift of $\ell$-torsion (beware that it is not unique: if $\tilde{P}$ is a cubical point of $\ell$-torsion, so is $\zeta\tilde{P}$ for $\zeta \in \mu_\ell$. Since $\lambda_0 = \lambda_1'$ in that case, $m\tilde{P}$ is also a canonical cubical point of $\ell$-torsion). Using the formula for $(u\ell + v)\tilde{P}$ above, we see that these two equations can be rewritten as a single equation $(\ell' + 1)\tilde{P} = -\ell'\tilde{P}$. Finding a canonical point of $\ell$-torsion thus corresponds to solving the equation $x^\ell = \beta$, and since $\beta^2 = e_{T,\ell}(P, P)$ and $\ell$ is odd, we see that we have an obstruction to finding a rational solution given by the self Tate pairing. Note however that even if a canonical point $\tilde{P}$ is given over an extension, $\tilde{P}^{\otimes \ell}$ is rational since $\beta$ is rational; and one can check (using the analytic
formulas of [LR22a]) that this point corresponds exactly to the unique element $g_P \in G(\ell D)$ which is symmetric of order $\ell$ (this is coherent with Example 2.7). By the discussion above, the other elements $g_P$ of order $\ell$ are given by the $\tilde{P}^{\otimes \ell}$ where $\tilde{P}$ satisfies the relaxed equation $(\ell' + 1)\tilde{P} = \zeta - \ell'\tilde{P}$ for some $\zeta \in \mu_\ell$.

Now assume that $\ell = 2m$ is even. The situation is more tricky. First, consider $\ell = 2$ and $P \in A[2]$; then $\lambda_0^2 = \lambda_1^{\prime 2}$, and we have two cases: $\lambda_0 = \lambda_1^{\prime}$ or $\lambda_0 = -\lambda_1^{\prime}$. One can check that the sign is given by $e_{D,\star}(P) = (-1)^{m_D(P)-m_D(0)}$ where $m_D(P)$ is the multiplicity of $D$ at $P$ (so for an elliptic curve and $D = (0_E)$, $e_{D,\star}(P) = 1$ for all $P \in E[2]$). In the later case, when $e_{D,\star}(P) = -1$, there are no canonical cubical points of 2-torsion above $P$. But when $e_{D,\star}(P) = 1$, the equation $2\tilde{P} = \tilde{0}$ has four solutions (over $\bar{k}$), and is enough for $\tilde{P}$ to be a canonical point of 2-torsion.

Now for the general even case $\ell = 2m = 2\ell' + 2$, $P \in A[\ell]$. Let $P_0 = mP$. We have $\lambda_0^2 = \lambda_1^{\ell'}$, so $\lambda_1^{\prime m} = \pm \lambda_0$. But $\lambda_0(P_0) = \lambda_0(P)$, $\lambda_1'(P_0) = \lambda_1'(P)^m$, so $\lambda_1^{\prime m} = e_{D,\star}(P_0)\lambda_0$. If $e_{D,\star}(P_0) = -1$, there are no canonical points of $\ell$-torsion. If $e_{D,\star}(P_0) = 1$, there are $2\ell$ canonical points of $\ell$-torsion, which can be found by solving the two equations $\ell'\tilde{P} = \tilde{0}$, $(\ell + 1)\tilde{P} = \tilde{P}$, or the single equation $(\ell' + 2)\tilde{P} = -\ell'\tilde{P}$. The solutions are parametrized by $\gamma^{2\ell} = \lambda_1^{\prime}$. These $2\ell$ solutions for $\mathcal{L}$ induces 2 distinct elements $\tilde{P}^{\otimes \ell}$, which correspond to the two symmetric elements above $P$ (necessarily of order $\ell$) in $G(\ell D)$. These two elements are rational if and only if $\lambda_1' = e_{T,\ell}(P, P)$ is a square; this is coherent with Example 2.7.

### 6.2.3. Radical isogenies.
We remark that the cubical arithmetic also extends to isogenies. Let $f: A \to B$ be an isogeny of abelian varieties, $M$ a line bundle on $B$ and $\mathcal{L}$ a line bundle on $A$ such that $f$ is a $\ell$-isogeny, i.e., $f^*M \simeq \mathcal{L}^\ell$. We fix such an isomorphism once and for all. We remark that a rigidification $\tilde{P}$ of $f^*M$ (hence via the isomorphism above, a rigidification of $\mathcal{L}^\ell$) at a point $P \in A$ corresponds to a rigidification of $M$ at the point $f(P)$, which we denote by $f(\tilde{P})$ (or rather $\tilde{f}(\tilde{P})$ where $\tilde{f}$ depends on the choice of isomorphism $f^*M \simeq \mathcal{L}^\ell$). It is convenient to normalise the neutral points so that $\tilde{f}(\tilde{0}_A) = \tilde{0}_B$.

Cubical points thus also give the correct conceptual framework to build algorithms for radical isogenies. It is known (see [CDV20; CDHV22] for dimension 1 and [CD21; LR22b] for the general case) that radical isogenies can be described by choices of $\ell$-th root of unity on Tate pairings.

Let us focus to the case of radical elliptic curve isogenies for simplicity. We start with $P \in E$ a point of $\ell$-torsion, and the goal is to find formulas not only for the isogeny $\varphi: E \to E' = E/\langle P \rangle$ but also to find a new point $P' \in E'$ of $\ell$-torsion such that the associated isogeny is not backtracking. It is not hard to see that such choices of $P'$ are in bijection with $\tilde{\varphi}^{-1}(P)$, and the theory of radical isogenies state that this fiber is in bijection with choices of $\ell$-th root of the non reduced Tate self pairing $e_{T,\ell}(P, P)$. We refer to [LR22b] for a geometric description of this isomorphism. Since cubical point arithmetic can be used to compute the Tate pairing, it is not surprising that they can explain radical isogenies, or more generally fibers of isogenies.

In brief, let $f: E_1 \to E_2$ be a cyclic isogeny of degree $\ell$, and $\tilde{f}$ the contragredient isogeny. Fixing a point $P \in \text{Ker } \tilde{f}$ gives through the Weil-Cartier pairing $e_f$ an isomorphism $\text{Ker } f \simeq \mu_\ell$, and if $Q \in E_2(Q)$, the (non reduced) Tate pairing $e_{T,f}(P, \tilde{Q}) = e_{T,\ell}(P, Q)$, seen as a $\mu_{\ell}$-étale torsor $\zeta^\ell = e_{T,\ell}(P, Q)$, is isomorphic to the fiber $f^{-1}(Q)$ which is a $\text{Ker } f$-torsor. (In other words, the Galoisian structure of $f^{-1}(Q)$ can be derived from the Galoisian structure of $\zeta^\ell = e_{T,\ell}(P, Q)$.) For an $\ell$-isogeny $f: A_1 \to A_2$ of principally polarised abelian varieties of dimension $g$, fixing a basis $P_1, \dots, P_g$ of $\text{Ker } \tilde{f}$ splits $\text{Ker } f$ as a product $\mu_\ell^g$, hence splits the
torsor $f^{-1}(Q)$ as a product of $\mu_\ell$-torsors, whose isomorphism classes are given by the torsors
$\zeta^\ell = e_{T,\ell}(P_i, Q)$.

Let us go back to an $\ell$-isogeny $f : E_1 \to E_2$ of elliptic curves for simplicity, and assume
$\ell$-odd. The symmetric divisor $\ell(0_{E_1})$ on $E_1$ descends (up to linear equivalence) to the
symmetric divisor $0_{E_2}$ on $E_2$, so by Mumford's theory corresponds to a lift $\tilde{K}$ of $\text{Ker} f$ to the
theta group $G(\ell(0_{E_1}))$ (in practice: this is the unique lift given by symmetric elements of
order $\ell$). Now let $P \in \text{Ker} \tilde{f}$ and $P' \in E_1[\ell]$ be a preimage of $P$ by $f$. By the discussion above,
a choice of cubical point $\tilde{P}$ (of level 1) for $P$ gives a choice of cubical point $\tilde{P}'$ (of level $\ell$) for
$P'$, but since $P'$ is in $E_1[\ell]$, the local rigidification of the line bundle associated to $\ell(0_E)$ at
$P'$ gives a global rigidification, hence an element $g_{P'} \in G(\ell(0_{E_1}))$, and conversely. So in the
other direction, since we have the canonical level subgroup $\tilde{K} \subset G(\ell(0_{E_1}))$ above $K = \text{Ker} f$,
we have a canonical cubical point of level $\ell$, for $P'$, hence a canonical cubical point $\tilde{P}$ of level 1
for $P$. One can check that this $\tilde{P}$ is symmetric and of order $\ell$ for the cubical arithmetic (in the
terminology above it is a canonical cubical point above $\tilde{P}$), and that replacing $P'$ by $P' + T$
with $T \in \text{Ker} f$ replace $\tilde{P}$ by $\zeta \tilde{P}$ for $\zeta = e_f(T, P)$. In particular, the $\ell$ canonical cubical points
$\tilde{P}$ above $P$ all come from a choice of $P' \in f^{-1}(P)$.

On the other hand, to compute a canonical cubical point $\tilde{P}$ above $P$, one can start with
an arbitrary cubical point $\tilde{P}$, write $\ell = 2\ell' + 1$, and compute the monodromy $\lambda_P$ between
$(\ell' + 1)\tilde{P}$ and $-\ell'\tilde{P}$. By Remark 3.13 and Section 4.2.5, this monodromy $\lambda_P$ is exactly the
square root of the non reduced Tate pairing $e_{T,\ell}(P, P)$. Canonical cubical points are then
of the form $\zeta \tilde{P}$, where $\zeta^{-\ell} = \lambda_P$. Since we have seen that we also had a bijection with the
preimages of $P$, this gives an explicit bijection between $\zeta^{-\ell} = \lambda_P = e_{T,\ell}(P, P)^{1/2}$ and
$f^{-1}(P)$.

To study the general fiber $f^{-1}(Q)$ for some point $Q \in E_2(k)$, we can proceed similarly.
Fix a $Q'$ in the preimage, and an arbitrary cubical point of level 1 for $Q$ (hence of level $\ell$ for
$Q'$). We have seen that each preimage of $P'$ of $P$ gives a level $\ell$ theta group element $g_{P'}$, hence
an action $g_{P'} \cdot \tilde{Q}'$, hence by descent a cubical point $P + \tilde{Q} = g_{P'} \cdot \tilde{Q}$ above $P + \tilde{Q}$. One can
check that $P + \tilde{Q}$ satisfy the equation $\ell P + \tilde{Q} = \tilde{Q}$, and that changing $P'$ by $P' + T$ changes
$P + \tilde{Q}$ to $e_f(T, P)P + \tilde{Q}$. Hence we get all $\ell$-cubical points $P + \tilde{Q}$ satisfying this equation.

On the other hand, we can also compute these cubical points on $E_2$ directly: take an
arbitrary $P + \tilde{Q}$, an arbitrary canonical cubical point $\tilde{P}$ so that $\ell \tilde{P} = \tilde{0}$, and compute the
monodromy $\ell \tilde{P} + \tilde{Q} = \lambda_P \tilde{Q}$. By the monodromy interpretation of the Tate pairing, we have
$\lambda_P = e_{T,\ell}(P, Q)$. The points $P + \tilde{Q}$ we search for are the ones given by $\zeta \tilde{P} + \tilde{Q}$ where $\zeta^{-\ell} =
\lambda_P$. Hence in this case, we have an explicit bijection between couples $(\tilde{Q}, \tilde{Q} + P)$ satisfying
$\tilde{Q} + \ell P = \tilde{Q}$; roots of the equation $\zeta^{-\ell} = \lambda_P = e_{T,\ell}(P, Q)$; and couples $(Q', Q' + P' + T)$
where $f(Q') = Q, f(P') = P$ are fixed and $T$ goes through $\text{Ker} f$.

The same strategy works for higher dimension. In practice, this allows to obtain explicit
formulas<fnref n="3" /> for fibers of isogenies between abelian varieties in the theta model, and also to
obtain explicit multiradical isogeny formulas in the theta model. The formulas can be derived
from the study of fibers of modular correspondances of [FLR11; Rob21a, § 5.2.2], then
using [LR10; LR15] to reinterpret the constants appearing in these formulas as suitable Tate
pairings, and then (a not yet published) adaptation of the descent level formula from [LR22a]
to go down in level without multiplying by $\ell$.

<footnote n="3">³These formulas, resulting from joint work with David Lubicz, were already promised more than one year ago in [Rob23c, Remark 5.19], but see footnote ²...</footnote>
**6.3. Supersingularity testing as a self pairing test.** Let $E/\mathbb{F}_q$ be an elliptic curve (with $j$-invariant different from 0, 1728), $q = p^2$. We would like to test if it is supersingular. In [Dol18], Doliskani proposes the following probabilistic supersingularity test.

Recall that on an elliptic curve with a Weierstrass equation, there are canonical polynomials $\phi_n, \omega_n, \psi_n$ (depending on $E$) such that if $P = (x_P, y_P)$, $nP = (\frac{\phi_n(x_P)}{\psi_n^2(x_P)}, \frac{\omega_n(x_P, y_P)}{\psi_n^3(x_P, y_P)})$. The polynomial $\psi_n$ is the $n$-th division polynomial: $\psi_n(x_P, y_P) = 0$ if and only if $P$ is a non trivial point of $n$-torsion. (If $n$ is odd, $\psi_n$ depends only on $x$).

Furthermore, there is a linear recurrence relation for these polynomials, which allows to compute $\psi_n(x_P, y_P)$ efficiently.

If the characteristic is $p$, then $\psi_p(x) = \tilde{\psi}_p(x)^p$, where $\tilde{\psi}_p(x)$ is a polynomial of degree $(p-1)/2$ if $E$ is ordinary, or $\pm 1$ if $E$ is supersingular. In particular, $E$ is supersingular if and only if $\psi_p^2 = 1$ [Dol18, Lemma 4; BGS22, Lemma 3].

Doliskani's supersingularity test is then to sample a random $x \in \mathbb{F}_q$, and test if $\psi_p(x)^2 = 1$. If it is not the case, we know that $E$ is ordinary. Otherwise, by the Schwartz-Zippel lemma [Dol18, Lemma 2], we know that $E$ is supersingular with high probability [Dol18, § 3]. In [Dol18, § 4] this test is further refined.

We now reinterpret this test as a pairing test. First, the division polynomials $\psi_n$ give precisely the elliptic nets of rank 1 (i.e., elliptic divisibility sequences) [Stao08, Theorem 1.2.1], and elliptic nets give pairings by [Stao08, § 17, § 18]. Here we will rather use the cubical point of view (this is essentially the same thing by Section 4.9.5).

Analytically we have $\Psi_n(z) = \frac{\sigma(nz)}{\sigma(z)^n}$. It follows from Sections 4.8 and 4.9 that if $\tilde{P} = (x(P), y(P), 1)$ is a level 3 cubical point given by level 3 affine Weierstrass coordinates, then $n\tilde{P} = ((\phi_n(x_P)\psi_n(x_P, y_P), \omega_n(x_P, y_P), \psi_n^3(x_P, y_P)))$. Likewise, if $\tilde{P}$ is a level 1 cubical point normalised by $Z_1(\tilde{P}) = 1$, then $Z_1(n\tilde{P}) = \psi_n(x_P, y_P)$. (As an aside, by Remark 4.24, using division polynomial for cubical arithmetic means that $\tilde{O}$ is normalised with respect to $-x/y$ rather than with respect to $x/y$).

So Doliskani's supersingularity test can be reinterpreted as follows: given $\tilde{P}$ normalised as above, is $Z_2(p\tilde{P}) := Z_1^2(p\tilde{P}) = 1$?

Now if $E$ is supersingular over $\mathbb{F}_q$ (with maximal endomorphism ring), then $E(\mathbb{F}_q) = (\mathbb{Z}/(p \pm 1)\mathbb{Z})^2$. We first sample a random point $x(P)$ on the Kummer line; $P$ lies either in $E$ or its quadratic twist, so $P$ is of order $p+1$ or $p-1$. In both cases, $x(pP) = x(P)$, this gives us a first test for a point order.

We can refine this test as follows: take $\tilde{P}$ a cubical point of level 1 normalised by $Z_1(\tilde{P}) = 1$, and of order $p \pm 1$, which we have checked via the equation $x(pP) = x(P)$. We want to compute the self Tate pairing $e_{p\pm 1}(P, P)$. By Theorem 4.19 and Remark 3.13, if $P$ is of order $p-1$, we can compute $e_p(P, P)$ by comparing $Z_1(p\tilde{P})$ and $Z_1(\tilde{P}) = 1$. In particular, $Z_1(p\tilde{P}) = 1$ (i.e. $\psi_p(P) = 1$) if and only if $e_{p-1}(P, P) = 1$. And if $P$ is of order $p+1$, we can compute $e_p(P, P)$ by comparing $Z_1(p\tilde{P})$ and $Z_1(-\tilde{P}) = -1$. In particular, $Z_1(p\tilde{P}) = -1$ (i.e. $\psi_p(P) = -1$) if and only if $e_{p+1}(P, P) = 1$. (We remark that if $\lambda \in \mathbb{F}_q = \mathbb{F}_{p^2}$, $p(\lambda \cdot \tilde{P}) = \lambda p^2 \cdot p\tilde{P} = p\tilde{P}$ so the Tate pairing is already reduced here.) So we can reframe Doliskani's supersingularity test as a self pairing test: is $e_{p\pm 1}(P, P) = 1$ for a randomly sampled $P$?

**Lemma 6.1.** Let $P \in E(\mathbb{F}_q)$ be of order $p \pm 1$. Then the reduced Tate self pairing $e_{T, p \pm 1}(P, P) = 1$ if and only if the elliptic curve quotient $E' = E/\langle P \rangle$ has its full $p \pm 1$-torsion rational over $\mathbb{F}_q$.
*Proof.* This follows from the geometric interpretation of the Tate pairing, see [Rob23c, Example 5.14]. □

In particular, if $E$ is supersingular, then since $E/\langle P \rangle$ has the same Galois structure than $E$ we always have $e_{T,p\pm 1}(P) = 1$: self pairings are trivial on a supersingular curve. (This can also be seen more directly as follows: if $P$ is of order $\ell$, $e_{T,\ell}(P) = e_{W,\ell}(P, \pi P' - P')$ where $\ell P' = P$, but $\pi$ is a scalar for a supersingular curve and the Weil pairing is alternate, hence the self pairing is trivial.) Lemma 6.1 thus gives an alternative proof that $\psi_p(x(P)) = \pm 1$ for a supersingular curve.

We can also use Lemma 6.1 and our reinterpretation of Doliskani's supersingularity test as a self pairing test to give a precise description of the points $P \in E(\mathbb{F}_q)$ such that Doliskani's test fails for an ordinary elliptic curve $E: P$ has to be of $p \pm 1$-torsion and the isogeneous curve $E/\langle P \rangle$ has to have fully rational $p \pm 1$-torsion. Using the group structure of isogeny volcanoes of ordinary curves, this allows to refine the probability of failure (depending on where $E$ is in the volcano).

In practice, in isogeny based cryptography we work with the Montgomery model. We can thus use the fast cubical ladder formulas from Section 5.2 to do our self pairing test: sample a random point $x(P)$, start with $\tilde{P}$ such that $X_2(\tilde{P}) = x(P)$, $Z_2(\tilde{P}) = 1$ in level 2 affine cubical coordinates, and compute whether $p\tilde{P} = \tilde{P}$. This test both that $P$ is of $p \pm 1$ torsion and that the self pairing is trivial (or rather its square, since we are using level 2 coordinates).

We recover precisely the fast supersingular test from [BGS22]. In that paper, the authors apply the above strategy, but using the usual projective Montgomery ladder to compute $p\tilde{P}$. As we have seen in Section 5.2, the standard projective ladder almost correctly computes the cubical arithmetic, and it is easy to keep track of a correcting factor to apply afterwards to obtain the correct cubical point $p\tilde{P}$. Now comparing Algorithms 5.4 and 5.5 with [BGS22, Algorithm 1], we can check that this correcting factor is precisely the factor from [BGS22, Proposition 2]. Since the authors of [BGS22] use this correcting factor in their fast supersingularity test [BGS22, Algorithm 6], they are really computing the true cubical exponentiation [BGS22]! (We remark that computing $1/x(P)$ and then directly using the formulas Algorithms 5.4 and 5.5 for the cubical ladder would save $1M$ by bit.)

**6.4. Monodromy leak and the DLP.** Our last non pairing application, and perhaps the most important one, is a new devastating side-channel attack against the Montgomery ladder on Montgomery curves, which we call the monodromy leak, a type of projective coordinate leak.

Let $P \in E(\mathbb{F}_q)$ be a point of $\ell$-torsion, with $\ell$ prime for simplicity. It is known since [NSS04] that a projective coordinate leak, that is leaking the individual projective coordinates $X(mP), Y(mP), Z(mP)$ which are computed during the scalar multiplication, where $mP = (X(mP) : Y(mP) : Z(mP))$, would yield information on the secret $m$. The attack of [NSS04] can only recover a few bits of $m$, so the attack was only used for attacking the signature scheme ECDSA. For ECDSA, obtaining a few bit of leakage for each signature, combined with lattices methods, allows to recover the full secret key. This attack was revisited in [AGB20] where the authors found that many implementation were still vulnerable to projective coordinate leaks, and extended it to the Montgomery ladder.

By contrast, our monodromy leak attack on the Montgomery ladder is much more devastating, since it allows to recover the full key via only one leak, by reduction to some DLPs in $\mathbb{F}_q^*.$
The main idea is as follows: first we can assume that $\mu_\ell \not\subset \mathbb{F}_q$, otherwise pairings already give an efficient reduction from DLPs to $E(\mathbb{F}_q)$ to DLPs in $\mathbb{F}_q^*$. We can also assume $\ell$ odd for simplicity. Under this hypothesis, for a point $P \in E[\ell](\mathbb{F}_q)$ there is only one rational cubical point $\tilde{P}$ which is still of $\ell$-torsion. We will call $\tilde{P}$ the canonical lift of $P$, and denote it also by $\hat{P}$. We can efficiently compute $\hat{P}$ by taking an arbitrary rational cubical lift $\tilde{P}$ and then computing $\hat{P} = u\tilde{P}$ where $u \equiv 1 \bmod \ell$ and $u \equiv 0 \bmod q-1$; such a $u$ exists because $\ell$ is prime to $q-1$ by assumption.

**Lemma 6.2.** If $\mu_\ell(\mathbb{F}_q) = 1$ (i.e. $\ell$ is coprime to $q-1$), there is a unique rational canonical cubical point of $\ell$-torsion $\hat{P}$. It suffices to check that $\ell\hat{P} = \tilde{0}$, and for any rational cubical point $\tilde{P}$ above $P$, $\hat{P} = u\tilde{P}$, where $u \equiv 1 \bmod \ell$ and $u \equiv 0 \bmod q-1$.

Furthermore, for any $m$, we have $m\hat{P} = m\hat{P}$

*Proof.* We defined canonical cubical points of $\ell$-torsion in Section 6.2; they satisfy $\ell\tilde{P} = \tilde{0}$, $(\ell+1)\tilde{P} = \tilde{P}$. Using the notations of Section 6.2.2, if $\ell\tilde{P} = \tilde{0}$ then $\lambda_0 = \beta^\ell = 1$, hence $\beta = 1$ since $\beta$ is rational and $\mu_\ell(\mathbb{F}_q) = 1$, so $\lambda'_1 = 1$, and in this case the equation $\ell\tilde{P}$ is sufficient to define a canonical point of $\ell$-torsion.

For a random cubical point $\tilde{P}$, if $\ell\tilde{P} = \lambda_0\tilde{0}$, then $\gamma\tilde{P}$ is a canonical point of $\ell$-torsion for $\gamma^\ell = \lambda_0^{-1}$. Since $\mu_\ell(\mathbb{F}_q) = 1$, $x \mapsto x^\ell$ is bijective in $\mathbb{F}_q^*$, so this equation has a unique solution in $\mathbb{F}_q^*$.

We have $u\hat{P} = \hat{P}$ because $u \equiv 1 \bmod \ell$ and $\hat{P}$ is a canonical cubical point of $\ell$-torsion. On the other hand, starting with an arbitrary rational cubical point $\tilde{P} = \lambda\hat{P}$, $\lambda \in \mathbb{F}_q$, we have $u\tilde{P} = \lambda^{u^2}u\hat{P}$, and since $u \equiv 0 \bmod q-1$, $\lambda^{u^2} = 1$. So $u\tilde{P} = \hat{P}$.

If $Q = mP$, taking an arbitrary rational $\tilde{P}$ and letting $\tilde{Q} = m\tilde{P}$, we obtain $\tilde{Q} = u\tilde{Q} = um\tilde{P} = mu\tilde{P} = m\hat{P}$. $\square$

(We remark that [Stao8, Chapter 19; LSo8] also compute an elliptic net representation of the level 1 canonical lift $\hat{P}$ for elliptic curves. Their formula to compute $\hat{P}$ is slightly different than the one above, but do give the same canonical lift.)

The relation $m\hat{P} = m\hat{P}$ can also be proved by unicity and the compatibility relations of the cubical arithmetic: in general, if $\tilde{P}$ is a canonical lift of $\ell$-torsion, so is $m\tilde{P}$ for any $m$. In our situation, the canonical lift is unique (if we want it to be rational), so by Corollary 4.33, we even have that if $[P_1, P_2, P_3, P_4; Q_1, Q_2, Q_3, Q_4]$ are points of $\ell$-torsion in Riemann position, so are their canonical lift of $\ell$-torsion.

Now the cubical arithmetic is a mix of elliptic curve arithmetic and $\mathbb{F}_q^*$ arithmetic, and we can use $\hat{P}$ to try to reduce to DLPs in $\mathbb{F}_q^*$. Namely, for the Montgomery ladder, one start with a normalised point $P = (x(P) : 1)$ in order to reduce the number of multiplications. Let us define $\tilde{P} = (x(P), 1)$, and assume we obtain a projective coordinate leak of $m.P = (X(mP), Z(mP))$. Now assume furthermore that $m.P$ was computed using the cubical formulas, in other words $X(mP), Z(mP)$ gives $m\tilde{P}$ exactly. Then we have $\tilde{P} = \lambda\hat{P}$ with $\lambda = 1/Z(\hat{P})$, and $m\tilde{P} = \lambda^{m^2}m\hat{P}$ by homogeneity.

Now we know $mP = (X(mP) : Z(mP))$ because it is public, we know $\tilde{P} = (x(P), 1)$ because $P$ is normalised, we also know $m\tilde{P} = (X(mP), Z(mP))$ because we have assumed we had a projective coordinate leak. We know the canonical lift $\hat{P}$, and we can also compute $m\hat{P}$ even if we don't know $m$, because the canonical lift $m\hat{P}$ of $mP$ is precisely equal to $m\hat{P}$.
So we obtain an equation of the form $\lambda^{m^2} = Z(\overline{m\hat{P}})/Z(\hat{m}\hat{P})$, where the only unknown is $m$. Solving a DLP in $\mathbb{F}_q^*$ followed by a square root allows to recover the possible values of $m$ modulo $q-1$, or more precisely modulo the order of $\lambda$.

In practice, we know that in the Montgomery ladder implementation, $m$ is as small as possible, so in particular is smaller than $\#E(\mathbb{F}_q) \le q+1+2\sqrt{q}$. We first solve the DLP in $\mathbb{F}_q^*$ to recover $m^2$ modulo the order $N$ of $\lambda$ (so that $N \mid q-1$). We then test all possible square roots, there are at most $2^t$ such square roots where $t$ is the number of distinct prime factors of $q-1$. Unless $q-1$ is very smooth, computing all square roots is not a bottleneck compared to the DLP in $\mathbb{F}_q^*$.

Now, for each possible square root $m$ modulo $N$, we need to test all possibilities $m + aN$ up to the upper bound $m \le \#E(\mathbb{F}_q)$. The number of possibilities is in $O(q/N)$, the smaller the order of $\lambda$ the more tries we need, but on the other hand the lower the probability to stumble on such a $\lambda$. By Merten's theorem, the average number of tries for a random $\lambda$ is in $O(\log q)$.

Note that the above method applies even if $\tilde{P}$ is not normalised, as long as we know both $\tilde{P}$ (say because the implementation is public) and $m\tilde{P}$ (because of a projective coordinate leak).

However, this only applies because we know that $m \le l$. If we don't know the projective coordinate leak $m\tilde{P}$, we can still take an arbitrary cubical point $\overline{m\tilde{P}}$ (well not quite arbitrary, we want to choose $Z(\overline{m\tilde{P}})$ such that $Z(\overline{m\tilde{P}})/Z(\hat{m}\hat{P})$ is in the group generated by $\lambda$). Then we have $\overline{m\tilde{P}} = M\tilde{P}$ for some $M \equiv m \bmod l$, but this time we only have the bound $M \le l(q-1)$. The cubical arithmetic allows to recover the value of $M$ modulo $q-1$ (or rather modulo $N$ the order of $\lambda$), which gives zero information on the value of $M$ modulo $l$ since $l$ is prime to $q-1$. In other words, the monodromy attack only applies if we know the number (or we can pin this number in a small interval) of loop around $l$ we did when computing $M.\tilde{P}$.

**Remark 6.3** (Biextension monodromy leak). We can also do a monodromy attack using biextension arithmetic, this requires more information but bypass the square root step. Namely, still with our hypothesis that $l$ is prime to $q-1$, given another rational point $Q$, there is a unique rational biextension element $\hat{g}_{P,Q}$ above $P$ of order $l$ for $\star_1$. We have $\hat{g}_{P,P}^{\star_1,m} = \hat{g}_{mP,P}$.

It follows that if through some projective coordinate leak we are able to recover both some biextension element $g_{P,P} = \lambda\hat{g}_{P,P}$ and $g_{P,P}^{\star_1,m} = \lambda^m\hat{g}_{mP,P}$, we can recover $m$ modulo the order $N \mid q-1$ of $\lambda$ by solving an equation $\lambda^m = C$.

In practice, this could happen if we have a projective coordinate leak of both $m\tilde{P}$, $(m+1)\tilde{P}$ (this is not unreasonable since the Montgomery ladder computes $mP$, $(m+1)P$).

Namely, we know that on input the ladder start with the normalised point $\tilde{P} = (x(P), 1)$; we compute $2\tilde{P}$ by a cubical doubling, and represent $g_{P,P}$ by the cubical points $[\tilde{P}, \tilde{P}; \tilde{0}, 2\tilde{P}]$. Then $g_{P,P}^{\star_1,m}$ is represented by $[\tilde{P}, m\tilde{P}; \tilde{0}, (m+1)\tilde{P}]$, so a projective leak of $m\tilde{P}$, $(m+1)\tilde{P}$ indeed gives us $g_{P,P}^{\star_1,m}$. (As an aside: $\hat{g}_{mP,P} = [\hat{P}, \hat{m}\hat{P}; \hat{0}, (\hat{m}+1)\hat{P}]$.)

Now, all the discussion above assumes that the Montgomery ladder is implemented using the cubical arithmetic; which is definitively *not* the case. But we saw in Section 1.2 that it is *very close* to the cubical ladder.
More precisely, the projective doubling formula is exactly the same as in Algorithm 5.4. However, the projective differential addition formula uses the equation

$$ 
\begin{align*}
U &= (X(P) - Z(P))(X(Q) + Z(Q)) \\
V &= (X(P) + Z(P))(X(Q) - Z(Q)) \\
X(P + Q) &= Z(P - Q)(U + V)^2 \\
Z(P + Q) &= X(P - Q)(U - V)^2
\end{align*}
$$

which differs from the cubical differential addition Algorithm 5.5 by the factor $4X(P - Q)Z(P - Q)$.

Taking this into account, we need to solve a slightly different degree two equations to recover $m$. Namely, let $[m\tilde{P}]$ be the cubical point computed by the usual Montgomery ladder. Then by [BGS22, Proposition 2], $[m\tilde{P}] = (4x(P))^{m2^{l(m)}-m}m\tilde{P}$ where $l(m)$ is the binary length of $m$ (see Section 6.3 for the link between the division polynomials in the statement of [BGS22, Proposition 2], and cubical arithmetic).

Using $\tilde{P} = \lambda_1\tilde{P}$, $Q = mP$, $m\tilde{P} = \lambda_1^{m^2}m\tilde{P}$, $[m\tilde{P}] = (4x(P))^{m2^{l(m)}-m}m\tilde{P}$, and since we know $P, Q$, hence $\tilde{P}, \tilde{Q}$, and also $\tilde{P}$ and $[m\tilde{P}]$ by assumption that we have a projective coordinate leak on the standard projective ladder, we can recover $\lambda_2$ such that $[m\tilde{P}] = \lambda_2 m\tilde{P}$ and we have the following equation where the only unknown is $m$:

$$ (22) \qquad (4x(P))^{m2^{l(m)}-m}\lambda_1^{m^2} = \lambda_2 $$

In practice, we also know the length $l(m)$ of $m$. So we fix $\zeta$ a primitive root of $\mathbb{F}_q^*$, and compute the dlp's with respect to $\zeta$: $\text{dlp}_\zeta(4x(P))$, $\text{dlp}_\zeta(\lambda_1)$, $\text{dlp}_\zeta(\lambda_2)$. Then by Equation (22), $m$ is a solution of the degree two equation:

$$ X^2(\text{dlp}_\zeta(\lambda_1) - \text{dlp}_\zeta(4x(P))) + 2^{l(m)} \text{dlp}_\zeta(4x(P))X - \text{dlp}_\zeta(\lambda_2) = 0. $$

We then proceed as before. The above discussion can thus be summarised in:

**Theorem 6.4.** Let $P = (X(P), Z(P))$ be a known public point of order $\ell$ on a Montgomery Kummer line associated to a Montgomery curve $E/\mathbb{F}_q$ (here we assume that we know not only $P$, but $X(P), Z(P)$, in practice $P$ is normalised via $Z(P) = 1$). Assume that $\ell$ is prime to $q-1$.

Let $m \le \ell$, and let $mP = (X(mP), Z(mP))$ as computed by the standard projective Montgomery ladder. Assume that we obtain a projective coordinate leak of $mP$, i.e., we not only know $x(mP) = X(mP)/Z(mP)$, but also $X(mP), Z(mP)$.

Let $u$ be the number of distinct prime factors of $q-1$. Let $\tilde{P} = (X(\tilde{P}), Z(\tilde{P}))$ be the unique canonical cubical rational point above $P$, $N | q-1$ be the multiplicative order of $Z(P)/Z(\tilde{P})$ and $v = (q-1)/N$.

Then one can recover $m$ by solving three discrete algorithms in $\mathbb{F}_q^*$ (two of which can be seen as a precomputation depending only on $P$, not $m$), followed by an algorithm polynomial in $\log q, 2^u$ and $v$.

Any constant time division algorithm to compute $x(mP) = X(mP)/Z(mP)$, such as the one employed by NaCL, prevents this attacks. Compared to [NSS04], we can recover $m$ fully from one leak, but we need to assume that we know $X(P), Z(P)$, not only $P$. So another protection is to mask $P$ by multiplying $X(P), Z(P)$ by some random scalar factor before doing the exponentiation.

For more details on the monodromy leak (and some wild speculations), we refer to [Rob23a, § 7.2] and to [Rob23b] for a toy implementation.
7. PERSPECTIVES

We have seen that biextensions provide a convenient framework to study pairings as monodromy information, and that the efficient representation of biextension elements provided by cubical points gave fast formulas to compute pairings on elliptic curves.

In this paper, we have mainly looked in Section 5 at cubical points on Montgomery Kummer lines represented by level 2 affine coordinates. It would be interesting to look at other models and also explore further a mix of projective and affine coordinates (notably in the Edwards model), as suggested in Remark 4.34. There is also potentially room for optimisation of the adding formulas in the double and add biextension algorithm for the Montgomery model. We have also seen in Section 6.1 that the cubical arithmetic could potentially be interesting for pairing based cryptography, not only for generic pairing computations. In fact, we saw that biextensions and cubical arithmetic shed new light (and gave new formulas) on Miller's algorithm.

In the other direction, one could ask if the refinements of Miller's algorithm from [BELL10; LL11; DZZZ22] (see Section 6.1.2) could also be applied in some sense to the cubical ladder. Notably so far in our explicit cubical arithmetic, we have only used doublings, differential additions and three way additions, but not the more general algebraic Riemann relations from Section 4.1. Although the general Riemann relations are induced by the squared torsor structure, and explicit formulas can be given via the differential addition formulas (see Section 4.9.3 for an example), expressing the general Riemann relations through differential addition formulas involve some divisions by 2, hence potentially points in field extensions, while using the formulas from Section 4.1 we can stay in the base field. Notably, the work of [BELL10] suggest looking at the following Riemann relations to compute the cubical exponentiation: [$(m_1+m_2)P+Q, -m_1P+Q, -m_2P+Q, -Q; -(m_1+m_2)P, m_1P, m_2P, 2Q$]. Maybe some other types of Riemann relations could also be useful.

We have focused in this article on computing pairings through monodromy (except in the setting of Remark 2.11. Furthermore, we work with biextension elements represented via cubical points, and we represent cubical points fully (or almost fully in level 2) through their cubical coordinates. This helps making the pairing formula complete, see Remarks 4.26 and 5.3: we can switch coordinates on the fly according to which point monodromy we want to compute. On the other hand, this means that for pairings, we compute redundant information: by Lemma 4.15 the cubical representation of a biextension element is not unique, and by Sections 4.5.2 and 4.6 cubical points are overkill to compute pairings, we just need cubical functions (a cubical point, represented by its cubical coordinates, encode a bunch of cubical functions at once). We already exploited the non unicity of the biextension representation in the double and add algorithm from Algorithm 4.1: the cubical points computed via this algorithm are not compatible with the cubical arithmetic, but they still give the correct biextension arithmetic. All in all this means that we have various ways to represent biextension elements (the evaluation representation, the cubical representation) and various ways to represent cubical points or cubical functions themselves, which makes our framework quite flexible, and there could still be potential for optimisations.

An interesting direction would also be to look at the arithmetic of a general biextension $X_f$ associated to an isogeny rather than just a polarisation, and see if it can help to compute the Weil-Cartier pairing $e_f$ and the generalised Tate pairing. It also could be helpful in deriving isogeny formulas; see Section 6.2 for why the cubical arithmetic can help to understand isogenies and their fibers.
One can also wonder whether this point of view can help about pairing inversion (we saw in Section 6.4 that it gave new insight on the dlp). From the point of view of biextensions, the reduced Tate pairing inversion is just the computation of a $q - 1$-th root in some theta group. But taking root of elements is not faster than a dlp for a generic group.

Going beyond pairings, we hope that Section 6 has convinced the reader that the notion of cubical points and cubical arithmetic is worthwhile to study, beyond just as a convenient way to work out the biextension arithmetic and computing the monodromy pairings. Indeed, some of our applications really need the full power of the cubical arithmetic. A cubical point is an amalgamation of an abelian variety point with some extra data coming from an action of $G_m$. Exploiting this action carefully through monodromy computations gave us back the usual pairings. One can wonder whether there could be some cryptographic applications which uses the cubical arithmetic directly, not just through pairing computations.

REFERENCES

[AGB20] A. C. Aldaya, C. P. García, and B. B. Brumley. "From A to Z: Projective coordinates leakage in the wild". In: *IACR Transactions on Cryptographic Hardware and Embedded Systems* (2020), pp. 428-453 (cit. on p. 76).

[BGS22] G. Banegas, V. Gilchrist, and B. Smith. "Efficient supersingularity testing over GF (p) and CSIDH key validation". In: *Mathematical Cryptology* 2.1 (2022), pp. 21-35 (cit. on pp. 67, 75, 76, 79).

[BBLP13] D. Bernstein, P. Birkner, T. Lange, and C. Peters. "ECM using Edwards curves". In: *Mathematics of Computation* 82.282 (2013), pp. 1139-1179 (cit. on p. 53).

[BDLS20] D. Bernstein, L. De Feo, A. Leroux, and B. Smith. "Faster computation of isogenies of large prime degree". In: *Algorithmic Number Theory Symposium (ANTS XIV)*. Vol. 4. 1. Mathematical Sciences Publishers, 2020, pp. 39-55. arXiv: 2003.10118. url: https://msp.org/obs/2020/4/p04.xhtml (cit. on p. 64).

[BBM79] P. Berthelot, L. Breen, and W. Messing. *Théorie de Dieudonné cristalline II*. Springer, 1979 (cit. on pp. 17, 18).

[BELL10] J. Boxall, N. El Mrabet, F. Laguillaumie, and D.-P. Le. "A variant of miller's formula and algorithm". In: *Pairing-Based Cryptography-Pairing 2010: 4th International Conference, Yamanaka Hot Spring, Japan, December 2010. Proceedings* 4. Springer. 2010, pp. 417-434 (cit. on pp. 1, 2, 9, 24, 67, 69, 70, 80).

[Bre83] L. Breen. *Fonctions thêta et théoreme du cube*. Vol. 980. Springer, 1983 (cit. on pp. 6, 9, 18, 19, 27, 35, 51, 54).

[BJ02] E. Brier and M. Joye. "Weierstraß elliptic curves and side-channel attacks". In: *International Workshop on Public Key Cryptography*. Springer. 2002, pp. 335-345 (cit. on p. 66).

[CLZ24] S. Cai, K. Lin, and C.-A. Zhao. "Pairing Optimizations for Isogeny-based Cryptosystems". In: *Cryptology ePrint Archive* (2024) (cit. on p. 9).

[CD21] W. Castryck and T. Decru. "Multiradical isogenies". In: *Cryptology ePrint Archive* (2021) (cit. on p. 73).

[CDHV22] W. Castryck, T. Decru, M. Houben, and F. Vercauteren. "Horizontal racewalking using radical isogenies". In: *International Conference on the Theory and Application of Cryptology and Information Security*. Springer. 2022, pp. 67-96 (cit. on p. 73).
[CDV20] W. Castryck, T. Decru, and F. Vercauteren. "Radical isogenies". In: *International Conference on the Theory and Application of Cryptology and Information Security*. Springer. 2020, pp. 493-519 (cit. on p. 73).

[CGF08] W. Castryck, S. Galbraith, and R. R. Farashahi. "Efficient arithmetic on elliptic curves using a mixed Edwards-Montgomery representation". In: *Cryptology ePrint Archive* (2008) (cit. on p. 3).

[CHM+23] W. Castryck, M. Houben, S.-P. Merz, M. Mula, S. v. Buuren, and F. Vercauteren. "Weak instances of class group action based cryptography via self-pairings". In: *Annual International Cryptology Conference*. Springer. 2023, pp. 762-792 (cit. on pp. 8, 9).

[Con04] B. Conrad. "Polarizations". 2004. URL: http://math.stanford.edu/~conrad/vigregroup/vigre04/polarization.pdf (cit. on pp. 8, 11).

[CR15] R. Cosset and D. Robert. "An algorithm for computing $(\ell, \ell)$-isogenies in polynomial time on Jacobians of hyperelliptic curves of genus 2". In: *Mathematics of Computation* 84.294 (Nov. 2015), pp. 1953-1975. DOI: 10.1090/S0025-5718-2014-02899-8. URL: http://www.normalesup.org/~robert/pro/publications/articles/niveau.pdf. HAL: hal-00578991, eprint: 2011/143. (Cit. on p. 71).

[CLN10] C. Costello, T. Lange, and M. Naehrig. "Faster pairing computations on curves with high-degree twists". In: *Public Key Cryptography–PKC 2010* (2010), pp. 224-242 (cit. on pp. 26, 27).

[DZZ23] Y. Dai, F. Zhang, and C.-a. Zhao. "Don't Forget Pairing-Friendly Curves with Odd Prime Embedding Degrees". In: *IACR Transactions on Cryptographic Hardware and Embedded Systems* 2023.4 (2023), pp. 393-419 (cit. on p. 9).

[DZZZ22] Y. Dai, Z. Zhou, F. Zhang, and C.-A. Zhao. "Software implementation of optimal pairings on elliptic curves with odd prime embedding degrees". In: *IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences* 105.5 (2022), pp. 858-870 (cit. on pp. 9, 24, 69, 80).

[Dol18] J. Doliskani. "On division polynomial PIT and supersingularity". In: *Applicable Algebra in Engineering, Communication and Computing* 29.5 (2018), pp. 393-407 (cit. on pp. 67, 75).

[FLR11] J.-C. Faugère, D. Lubicz, and D. Robert. "Computing modular correspondences for abelian varieties". In: *Journal of Algebra* 343.1 (Oct. 2011), pp. 248-277. DOI: 10.1016/j.jalgebra.2011.06.031. arXiv: 0910.4668 [cs.SC]. URL: http://www.normalesup.org/~robert/pro/publications/articles/modular.pdf. HAL: hal-00426338. (Cit. on p. 74).

[FK22] E. V. Flynn and K. Khuri-Makdisi. "An analog of the Edwards model for Jacobians of genus 2 curves". In: *arXiv preprint arXiv:2211.01450* (2022) (cit. on p. 53).

[GL08] S. Galbraith and X. Lin. "Computing Pairings Using x-Coordinates Only". In: *Designs, Codes and Cryptography* (2008). to appear (cit. on pp. 3, 50).

[GJMRV11] R. R. Goundar, M. Joye, A. Miyaji, M. Rivain, and A. Venelli. "Scalar multiplication on Weierstraß elliptic curves from Co-Z arithmetic". In: *Journal of cryptographic engineering* 1 (2011), pp. 161-176 (cit. on p. 66).

[Gro72] A. Grothendieck. *Groupes de Monodromie en Géométrie Algébrique: SGA 7*. Springer-Verlag, 1972 (cit. on pp. 4, 17-20).
[Ham20] M. Hamburg. "Faster Montgomery and double-add ladders for short Weierstrass curves". In: *IACR Transactions on Cryptographic Hardware and Embedded Systems* (2020), pp. 189-208 (cit. on p. 66).

[Koh11a] D. Kohel. "Addition law structure of elliptic curves". In: *Journal of Number Theory* 131 (2011), pp. 894-919 (cit. on p. 53).

[Koh11b] D. Kohel. "Arithmetic of split Kummer surfaces: Montgomery endomorphism of Edwards products". In: *Coding and Cryptology*. Springer, 2011, pp. 238-245 (cit. on p. 53).

[LS08] K. E. Lauter and K. E. Stange. "The elliptic curve discrete logarithm problem and equivalent hard problems for elliptic divisibility sequences". In: *International Workshop on Selected Areas in Cryptography*. Springer. 2008, pp. 309-327 (cit. on pp. 8, 77).

[LL11] D.-P. Le and C.-L. Liu. "Refinements of Miller's algorithm over Weierstrass curves revisited". In: *The Computer Journal* 54.10 (2011), pp. 1582-1591 (cit. on pp. 70, 80).

[LWXZ23] K. Lin, W. Wang, Z. Xu, and C.-A. Zhao. "A faster software implementation of sqisign". In: *Cryptology ePrint Archive* (2023) (cit. on p. 9).

[LR10] D. Lubicz and D. Robert. "Efficient pairing computation with theta functions". In: ed. by G. Hanrot, F. Morain, and E. Thomé. Vol. 6197. Lecture Notes in Comput. Sci. 9th International Symposium, Nancy, France, ANTS-IX, July 19-23, 2010, Proceedings. Springer-Verlag, July 2010. DOI: 10.1007/978-3-642-14518-6_21. URL: http://www.normalesup.org/~robert/pro/publications/articles/pairings.pdf. Slides: 2010-07-ANTS-Nancy.pdf (30min, International Algorithmic Number Theory Symposium (ANTS-IX), July 2010, Nancy), HAL: hal-00528944. (Cit. on pp. 6, 9, 51, 74).

[LR12] D. Lubicz and D. Robert. "Computing isogenies between abelian varieties". In: *Compositio Mathematica* 148.5 (Sept. 2012), pp. 1483-1515. DOI: 10.1112/S0010437X12000243. arXiv: 1001.2016 [math.AG]. URL: http://www.normalesup.org/~robert/pro/publications/articles/isogenies.pdf. HAL: hal-00446062. (Cit. on p. 71).

[LR15] D. Lubicz and D. Robert. "A generalisation of Miller's algorithm and applications to pairing computations on abelian varieties". In: *Journal of Symbolic Computation* 67 (Mar. 2015), pp. 68-92. DOI: 10.1016/j.ffa.2016.01.009. URL: http://www.normalesup.org/~robert/pro/publications/articles/optimal.pdf. HAL: hal-00806923, eprint: 2013/192. (Cit. on pp. 6, 9, 24, 50, 51, 74).

[LR16] D. Lubicz and D. Robert. "Arithmetic on Abelian and Kummer Varieties". In: *Finite Fields and Their Applications* 39 (May 2016), pp. 130-158. DOI: 10.1016/j.ffa.2016.01.009. URL: http://www.normalesup.org/~robert/pro/publications/articles/arithmetic.pdf. HAL: hal-01057467, eprint: 2014/493. (Cit. on pp. 3, 49, 51, 53, 66).

[LR22a] D. Lubicz and D. Robert. "Fast change of level and applications to isogenies". In: *Research in Number Theory (ANTS XV Conference)* 9.1 (Dec. 2022). DOI: 10.1007/s40993-022-00407-9. URL: http://www.normalesup.org/~robert/pro/publications/articles/change_level.pdf. HAL: hal-03738315. (Cit. on pp. 33, 71-74).

[LR22b] D. Lubicz and D. Robert. "Multiradical isogenies in the theta model". Sept. 2022. In preparation. (Cit. on p. 73).
[Mel07] N. Meloni. "New point addition formulae for ECC applications". In: *Arithmetic of Finite Fields: First International Workshop, WAIFI 2007, Madrid, Spain, June 21-22, 2007. Proceedings* 1. Springer. 2007, pp. 189-201 (cit. on p. 66).

[Mor85] L. Moret-Bailly. *Pinceaux de variétés abéliennes*. Société mathématique de France, 1985 (cit. on pp. 6, 9, 19, 27, 29, 31, 35, 52).

[Mum66] D. Mumford. "On the equations defining abelian varieties. I". In: *Invent. Math.* 1 (1966), pp. 287-354 (cit. on pp. 10, 12, 13, 29, 50, 52, 71, 72).

[Mum69] D. Mumford. "Bi-extensions of formal groups". In: *Algebraic geometry* 307-322 (1969) (cit. on pp. 4, 17).

[Mum70] D. Mumford. *Abelian varieties*. Tata Institute of Fundamental Research Studies in Mathematics, No. 5. Published for the Tata Institute of Fundamental Research, Bombay, 1970, pp. viii+242 (cit. on pp. 10, 11).

[Mum83] D. Mumford. *Tata lectures on theta I*. Vol. 28. Progress in Mathematics. With the assistance of C. Musili, M. Nori, E. Previato and M. Stillman. Boston, MA: Birkhäuser Boston Inc., 1983, pp. xiii+235. ISBN: 3-7643-3109-7 (cit. on p. 51).

[Mum84] D. Mumford. *Tata lectures on theta II*. Vol. 43. Progress in Mathematics. Jacobian theta functions and differential equations, With the collaboration of C. Musili, M. Nori, E. Previato, M. Stillman and H. Umemura. Boston, MA: Birkhäuser Boston Inc., 1984, pp. xiv+272. ISBN: 0-8176-3110-0 (cit. on p. 51).

[NSS04] D. Naccache, N. P. Smart, and J. Stern. "Projective coordinates leak". In: *Advances in Cryptology-EUROCRYPT 2004: International Conference on the Theory and Applications of Cryptographic Techniques, Interlaken, Switzerland, May 2-6, 2004. Proceedings* 23. Springer. 2004, pp. 257-267 (cit. on pp. 76, 79).

[OKUO11] N. Ogura, N. Kanayama, S. Uchiyama, and E. Okamoto. "Cryptographic pairings based on elliptic nets". In: *Advances in Information and Computer Security: 6th International Workshop, IWSEC 2011, Tokyo, Japan, November 8-10, 2011. Proceedings* 6. Springer. 2011, pp. 65-78 (cit. on pp. 54, 56).

[Rei23] K. Reijnders. "Effective Pairings in Isogeny-Based Cryptography". In: *International Conference on Cryptology and Information Security in Latin America*. Springer. 2023, pp. 109-128 (cit. on pp. 1, 9).

[Rob21a] D. Robert. "Efficient algorithms for abelian varieties and their moduli spaces". HDR thesis. Université Bordeaux, June 2021. URL: http://www.normalesup.org/~robert/pro/publications/academic/hdr.pdf. Slides: 2021-06-HDR-Bordeaux.pdf (1h, Bordeaux). (Cit. on pp. 3, 14, 23, 50, 70, 71, 74).

[Rob21b] D. Robert. *General theory of abelian varieties and their moduli spaces*. Mar. 2021. URL: http://www.normalesup.org/~robert/pro/publications/books/avtheory.pdf. Draft version. (Cit. on pp. 14, 50).

[Rob23a] D. Robert. "Improving the arithmetic of Kummer lines". Aug. 2023. URL: http://www.normalesup.org/~robert/pro/publications/notes/2023-11-kummer_lines.pdf (cit. on p. 79).

[Rob23b] D. Robert. "Kummer Line". Toolbox for computing on Kummer lines. Oct. 2023. URL: https://gitlab.inria.fr/roberdam/kummer-line (cit. on pp. 2, 9, 65, 69, 79).

[Rob23c] D. Robert. "The geometric interpretation of the Tate pairing and its applications". Feb. 2023. URL: http://www.normalesup.org/~robert/pro/
publications/articles/geometric_tate_pairing.pdf. eprint: 2023/177, HAL: hal-04295743v1. (Cit. on pp. 13, 14, 74, 76).

[RS24] D. Robert and N. Sarkis. "Computing 2-isogenies between Kummer lines". In: *IACR Communications in Cryptology* 1 (1 Jan. 2024). URL: http://www.normalesup.org/~robert/pro/publications/articles/kummer_isogenies.pdf. eprint: 2024/037. (Cit. on p. 13).

[Sil86] J. H. Silverman. *The arithmetic of elliptic curves*. Vol. 106. Graduate Texts in Mathematics. Corrected reprint of the 1986 original. New York: Springer-Verlag, 1986, pp. xii+400. ISBN: 0-387-96203-4 (cit. on p. 54).

[Sta08] K. Stange. "Elliptic nets and elliptic curves". PhD thesis. Brown University, 2008. URL: https://repository.library.brown.edu/studio/item/bdr:309/PDF/ (cit. on pp. 4–6, 18, 20, 32, 46, 48, 54, 56, 75, 77).

[Sta11] K. Stange. "Elliptic nets and elliptic curves". In: *Algebra & Number Theory* 5.2 (2011), pp. 197–229 (cit. on pp. 4, 6, 9).

[Tra14] C. Tran. "Formules d'addition sur les jacobiennes de courbes hyperelliptiques: application à la cryptographie". PhD thesis. Rennes 1, 2014. URL: https://www.theses.fr/185903150 (cit. on pp. 6, 56).

[Ver10] F. Vercauteren. "Optimal pairings". In: *IEEE Transactions on Information Theory* 56.1 (2010), pp. 455–461 (cit. on p. 25).

INRIA BORDEAUX-SUD-OUEST, 200 AVENUE DE LA VIEILLE TOUR, 33405 TALENCE CEDEX FRANCE

Email address: damien.robert@inria.fr

URL: http://www.normalesup.org/~robert/

INSTITUT DE MATHÉMATIQUES DE BORDEAUX, 351 COURS DE LA LIBERATION, 33405 TALENCE CEDEX FRANCE