In this section we compute, in detail, various forms of the four variable addition formulae for the quartic elliptic curve.
By quartic elliptic curve we mean the curve \begin{equation} \label{eq:quartic} y^2 = ax^4 + bx^3 + cx^2 + dx + e \end{equation} For convenience, let $R(x)$ denote the polynomial on the RHS of \eqref{eq:quartic}. As described in a previous section the four variable, level one, addition formula for the general quartic is given by \begin{equation} \label{eq:add4} \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x_4 & x_4^2 & y_4 \\ \end{vmatrix} \space = \space 0 \end{equation}
Our objective is to algebraically manipulate \eqref{eq:add4}, subject to the conditions $\{y_i^2 - R(x_i) = 0 : i=1,2,3,4\}$, into an addition formula of the form \begin{equation} \label{eq:XY} x_4 \space = \space X(x_1,x_2,x_3,y_1,y_2,y_3), \qquad y_4 \space = \space Y(x_1,x_2,x_3,y_1,y_2,y_3) \end{equation} where $X$ and $Y$ are rational functions.
Key Idea Part 1
The key idea in computing $x_4$ is that if you eliminate $y_4$ between the equations $y_4^2 = R(x_4)$ and \eqref{eq:add4} you get a polynomial of degree 4 in $x_4$. And three of roots of that polynomial $x_1,x_2,x_3$, are already known! Therefore we can infer the existance an identity of the form \begin{equation} \label{eq:key} \resultant_{y_4}\left(y_4^2 - R(x_4), \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x_4 & x_4^2 & y_4 \\ \end{vmatrix} \right) \space = \space (x_4 - x_1)(x_4 - x_2)(x_4 - x_3)(Dx_4 - N) \space + \space \sum_{i=1}^3 Q_i\big(y_i^2 - R(x_i)\big) \end{equation} where $D,N$ are polynomials in $x_1,x_2,x_3,y_1,y_2,y_3$ and $Q_1,Q_2,Q_3$ are polynomials in $x_1,x_2,x_3,x_4$ (because the resultant is of degree at most 2 in $y_1,y_2,y_3$). The value of $x_4$ is then given by \begin{equation} \label{eq:x4} x_4 \space = \space \frac N D \end{equation} We can confirm inference \eqref{eq:key} by computing the Normal Form of it's left hand side with respect to the Gröbner basis $\left\{ y_i^2 - R(x_i) : i=1,2,3 \right\}$ using CAS. Better still it is within reach of manual calculation ...
Manual Computation Of The Resultant Identity
Let $M_1,M_2,M_3,M_4$ be the cofactors of the last column of \eqref{eq:add4} \begin{equation} M_1 \space = \space - \begin{vmatrix} 1 & x_2 & x_2^2 \\ 1 & x_3 & x_3^2 \\ 1 & x_4 & x_4^2 \\ \end{vmatrix} \qquad M_2 \space = \space \begin{vmatrix} 1 & x_1 & x_1^2 \\ 1 & x_3 & x_3^2 \\ 1 & x_4 & x_4^2 \\ \end{vmatrix} \qquad M_3 \space = \space - \begin{vmatrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ 1 & x_4 & x_4^2 \\ \end{vmatrix} \qquad M_4 \space = \space \begin{vmatrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ 1 & x_3 & x_3^2 \\ \end{vmatrix} \end{equation} and $C_1,C_2,C_3,C_4$ be the omitted factors ie. \begin{equation} C_1 \space = \space (x_1-x_2)(x_1-x_3)(x_1-x_4) \qquad C_2 \space = \space (x_2-x_1)(x_2-x_3)(x_2-x_4) \qquad C_3 \space = \space (x_3-x_1)(x_3-x_2)(x_3-x_4) \qquad C_4 \space = \space (x_4-x_1)(x_4-x_2)(x_4-x_3) \end{equation} Then we have the algebraic identity
Utilising the identity $ \enspace M_1 + M_2 + M_3 + M_4 \space = \space 0 \enspace $ and then the identities $ \enspace M_1M_2 \space = \space C_3C_4, \enspace M_1M_3 \space = \space C_2C_4, \enspace M_2M_3 \space = \space C_1C_4 \enspace $ we can compute \begin{equation*} \begin{aligned} LHS \enspace &= \enspace - \left(M_1 y_1 + M_2 y_2 + M_3 y_3 + M_4 y_4\right)\left(M_1 y_1 + M_2 y_2 + M_3 y_3 - M_4 y_4\right) \\\\ &= \enspace M_1 (M_2+M_3+M_4) y_1^2 \enspace + \enspace M_2 (M_1+M_3+M_4) y_2^2 \enspace + \enspace M_3 (M_1+M_2+M_4) y_3^2 \enspace + \enspace M_4^2 y_4^2 \enspace - \enspace 2 M_1 M_2 y_1y_2 \enspace - \enspace 2 M_1 M_3 y_1 y_3 \enspace - \enspace 2 M_2 M_3 y_2 y_3 \\\\ &= \enspace M_4 \big( M_1 y_1^2 + M_2 y_2^2 + M_3 y_3^2 + M_4 y_4^2 \big) \enspace + \enspace M_1 M_2 (y_1-y_2)^2 \space + \space M_1 M_3 (y_1-y_3)^2 \space + \space M_2 M_3 (y_2-y_3)^2 \\\\ &= \enspace M_4 \sum_{i=1}^4 M_iR(x_i) \enspace + \enspace C_3 C_4 (y_1-y_2)^2 \space + \space C_2 C_4 (y_1-y_3)^2 \space + \space C_1 C_4 (y_2-y_3)^2 \enspace + \enspace M_4 \sum_{i=1}^4 M_i\big[y_i^2-R(x_i)\big] \end{aligned} \end{equation*} and we are almost there. Utilising the depleted Vandermonde determinant formula we have \begin{equation*} \sum_{i=1}^4 M_iR(x_i) \enspace = \enspace \begin{vmatrix} 1 & x_1 & x_1^2 & R(x_1) \\ 1 & x_2 & x_2^2 & R(x_2) \\ 1 & x_3 & x_3^2 & R(x_3) \\ 1 & x_4 & x_4^2 & R(x_4) \\ \end{vmatrix} \enspace = \enspace a \begin{vmatrix} 1 & x_1 & x_1^2 & x_1^4 \\ 1 & x_2 & x_2^2 & x_2^4 \\ 1 & x_3 & x_3^2 & x_3^4 \\ 1 & x_4 & x_4^2 & x_4^4 \\ \end{vmatrix} \enspace + \enspace b \begin{vmatrix} 1 & x_1 & x_1^2 & x_1^3 \\ 1 & x_2 & x_2^2 & x_2^3 \\ 1 & x_3 & x_3^2 & x_3^3 \\ 1 & x_4 & x_4^2 & x_4^3 \\ \end{vmatrix} \enspace = \enspace M_4 C_4 \big[a(x_1+x_2+x_3+x_4) \space + \space b\big] \end{equation*} which we can substitute in, then factor out $C_4$ to give \eqref{eq:add4identity}
Substituting $y_4 = \sqrt{R(x_4)}$ into \eqref{eq:add4identity}, negating and noting that the first term in square brackets is a polynomial of degree 1 in $x_4$, gives \eqref{eq:key} in explicit form.
Computation Of Formula For $x_4$
From \eqref{eq:add4identity} we have \begin{align*} N - D x_4 \space &= \space M_4^2 \big[a(x_1+x_2+x_3+x_4) + b\big] \space + \space \big[C_3 (y_1-y_2)^2 + C_2 (y_1-y_3)^2 + C_1 (y_2-y_3)^2\big] \\ &= \space (x_1-x_2)^2(x_1-x_3)^2(x_2-x_3)^2\big[a(x_1+x_2+x_3+x_4) + b\big] \space + \space (x_3-x_1)(x_3-x_2)(x_3-x_4)(y_1-y_2)^2 \\ &\qquad \qquad \qquad \qquad \qquad \qquad + \space (x_2-x_1)(x_2-x_3)(x_2-x_4)(y_1-y_3)^2 \space + \space (x_1-x_2)(x_1-x_3)(x_1-x_4)(y_2-y_3)^2 \end{align*} and can isolate $N$ and $D$, then $x_4 = N / D$ gives
\begin{equation} \label{eq:firstx4} x_4 \space = \space \frac {(x_1-x_2)^2(x_1-x_3)^2(x_2-x_3)^2\big[a(x_1+x_2+x_3) + b\big] \space + \space x_3(x_3-x_1)(x_3-x_2)(y_1-y_2)^2 \space + \space x_2(x_2-x_1)(x_2-x_3)(y_1-y_3)^2 \space + \space x_1(x_1-x_2)(x_1-x_3)(y_2-y_3)^2} {-a(x_1-x_2)^2(x_1-x_3)^2(x_2-x_3)^2 \space + \space (x_3-x_1)(x_3-x_2)(y_1-y_2)^2 \space + \space (x_2-x_1)(x_2-x_3)(y_1-y_3)^2 \space + \space (x_1-x_2)(x_1-x_3)(y_2-y_3)^2} \end{equation}
Addition Formula As Ratio Of Resultants
Let \begin{equation} \label{eq:G} G(x,y) \enspace = \enspace \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x & x^2 & y \\ \end{vmatrix} \end{equation} then from \eqref{eq:key} we have \begin{equation} \label{eq:residx4} \mathfrak{R}(x) \space = \space \resultant_y(F(x,y),G(x,y)) \space = \space D(x-x_1)(x-x_2)(x-x_3)(x-x_4) \end{equation} and since $D$ is independent of $x_4,y_4$ by evaluating \eqref{eq:residx4} at $x=0$ and $x=\infty$ we have
\begin{equation} \label{eq:resadd} x_4 \enspace = \enspace \frac 1 {x_1 x_2 x_3} \cdot \frac {\mathfrak{R}(0)} {\lim\limits_{t \rightarrow \infty} x^{-4} \thinspace \mathfrak{R}(x)} \qquad\qquad \end{equation}
Computation Of Addition Formula For $x_4$
To compute the resultant $\mathfrak{R}$, multiply together the determinants evaluated at the two $y_4$-conjugates for $x_4$, that is the $y$-roots of $y^2 - R(x_4) = 0$, giving using \eqref{eq:key} \begin{equation} \label{eq:X} \mathfrak{R} \space = \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x_4 & x_4^2 & y_4 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x_4 & x_4^2 & -y_4 \\ \end{vmatrix} \space = \space (x_4 - x_1)(x_4 - x_2)(x_4 - x_3)(Dx_4 - N) \end{equation} when all four points $(x_i,y_i)$ lie on the curve $y^2 = R(x)$.
Then obtain an expression for $N$ by evaluating \eqref{eq:X} at $(x_4,y_4) = (0,\sqrt{e})$ and an expression for $D$ by evaluating \eqref{eq:X} at $(x_4,y_4) = \big(t,\sqrt{R(t)}\big)$ as $t \rightarrow \infty$ to give \begin{equation} \label{eq:NDx4} x_1 x_2 x_3 N \space = \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & 0 & 0 & \sqrt{e} \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & 0 & 0 & -\sqrt{e} \\ \end{vmatrix}, \qquad D \space = \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 0 & 1 & \sqrt{a} \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 0 & 1 & -\sqrt{a} \\ \end{vmatrix} \end{equation}
Computation Of Addition Formula For $y_4$
Similarly to compute $y_4$, compute the resultant $\hat{\mathfrak{R}}$ by multiplying together the determinants evaluated at the four $x_4$-conjugates say $\chi_1,\chi_2,\chi_3,\chi_4$ of the $x$-roots of $R(x) - y_4^2 = 0$, giving \begin{equation} \label{eq:Y} \hat{\mathfrak{R}} \space = \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & \chi_1 & \chi_1^2 & y_4 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & \chi_2 & \chi_2^2 & y_4 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & \chi_3 & \chi_3^2 & y_4 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & \chi_4 & \chi_4^2 & y_4 \\ \end{vmatrix} \space = \space (y_4-y_1)(y_4-y_2)(y_4-y_3)(\hat{D}y_4 - \hat{N}) \end{equation} Let $e_1,e_2,e_3,e_4$ be the roots of $R(x)=0$. Obtain an expression for $\hat{N}$ by evaluating \eqref{eq:Y} at $(e_1,0)$ so that $\chi_i=e_i$. And obtain a formula for $\hat{D}$ by evaluating it at $\big(t,t^2\sqrt{a}\big)$ as $t \rightarrow \infty$, so that $\chi_i \rightarrow (-1)^{\frac{i-1} 2} t$, giving \begin{equation} \label{eq:NDy4} y_1 y_2 y_3 \hat{N} \space = \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & e_1 & e_1^2 & 0 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & e_2 & e_2^2 & 0 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & e_3 & e_3^2 & 0 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & e_4 & e_4^2 & 0 \\ \end{vmatrix},\qquad a^2 \hat{D} \space = \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 0 & 1 & \sqrt{a} \\ \end{vmatrix}^2 \thinspace \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 0 & 1 & -\sqrt{a} \\ \end{vmatrix}^2 \end{equation}
Four Variable Level 2 Addition Formula
From \eqref{eq:NDx4} and \eqref{eq:NDy4} we get an expression for $x_4$ and $y_4$ in terms of products of determinants
\begin{equation} \label{eq:x4y4} x_4 \space = \space \frac 1 {x_1 x_2 x_3} \frac {\begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & 0 & 0 & \sqrt{e} \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & 0 & 0 & -\sqrt{e} \\ \end{vmatrix}} {\begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 0 & 1 & \sqrt{a} \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 0 & 1 & -\sqrt{a} \\ \end{vmatrix}}, \qquad\qquad y_4 \space = \space \frac {a^2} {y_1 y_2 y_3} \frac {\begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & e_1 & e_1^2 & 0 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & e_2 & e_2^2 & 0 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & e_3 & e_3^2 & 0 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & e_4 & e_4^2 & 0 \\ \end{vmatrix}} {\begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 0 & 1 & \sqrt{a} \\ \end{vmatrix}^2 \thinspace \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 0 & 1 & -\sqrt{a} \\ \end{vmatrix}^2} \end{equation}
Observe that, due to the symmetries in \eqref{eq:x4y4} both $x_4$ and $y_4$ are rational functions of $x_1,x_2,x_3,y_1,y_2,y_3$ and the coefficients $a,b,c,d,e$ of the quartic curve.
RATIONAL FORMULA
CUBIC FORMULA
CIRCLE FORMULA
INTEGRAL FORMULA
ALTERNATE $x_4$ FORMULA
Alternate Addition Formula For $x_4$
You can get \eqref{eq:x4y4} by letting $x_5=0$ and $x_6 \rightarrow \infty$ in \begin{equation} \frac {(x_5-x_1)(x_5-x_2)(x_5-x_3)(x_5-x_4)} {(x_6-x_1)(x_6-x_2)(x_6-x_3)(x_6-x_4)} \space = \space \frac {\begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x_5 & x_5^2 & y_5 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x_5 & x_5^2 & -y_5 \\ \end{vmatrix}} {\begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x_6 & x_6^2 & y_6 \\ \end{vmatrix} \space \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x_6 & x_6^2 & -y_6 \\ \end{vmatrix}} \end{equation} By letting $x_5 \rightarrow x_1$ and $x_6 \rightarrow x_2$ you get \begin{equation} -\frac {(x_1-x_3)(x_1-x_4)} {(x_2-x_3)(x_2-x_4)} \space = \space \frac {\begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 2y_1 & 4x_1y_1 & R'(x_1) \\ \end{vmatrix}} {\begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 0 & 2y_2 & 4x_2y_2 & R'(x_2) \\ \end{vmatrix}} \end{equation} then what?
Geometric Interpretation of Addition Formula
As described in a previous section the three variable addition formula for the general cubic curve is given by \begin{equation} \label{eq:add3} \begin{vmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \\ \end{vmatrix} \space = \space 0 \end{equation} This addition formula for has a very simple geometric interpretation. Draw the cubic curve and a straight line on the cartesian plane $\R^2$. If the line intersects the curve at three points, then those three points satisfy the three variable addition formulae. This is a simple consequence of the fact that three points $(x_1,y_1) \space \ldots \space (x_3,y_3)$ lie on a straight line if and only if \eqref{eq:add3} holds.
This is often expressed in terms of finding rational points on the curve. If you have two rational points on a cubic curve, then to find a third draw a straight through the first two - the third point of intersection with the cubic curve will be another rational point. That does not work for a quartic curve, but ...
Four points $(x_1,y_1) \space \ldots \space (x_4,y_4)$ lie on a parabola $A + Bx + Cx^2 + Dy = 0$ if and only if \eqref{eq:add4} holds. So if you draw a parabola and it intersects a quartic curve at four points, then those points satisfy the four variable addition formulae. Therefore if you have three rational points on a quartic curve, you can find a fourth by drawing a parabola through the first three - the fourth point of intersection will be another rational point given by equation \eqref{eq:x4y4}
Rationalising the Denominators
Now for the "bad" news. Formulae \eqref{eq:x4y4} do not simply reduce to the usual classic three variable addition formulae when $(x_3,y_3)$ is set equal to an obvious rational point on the curve. For example when $(x_3,y_3)=(0,-1)$ on the curve $y^2 = 1 - x^4$ we do not immediately get the Euler's classic addition formula for the lemniscatic integral.
The reason for this is that in the classic formulae the denominators are "rationalised", that is they are polynomials in $x_1,x_2,x_3$ only. It is in principle easy to convert \eqref{eq:x4y4} into that form by multiplying the numerator and denominator by conjugates of the denominator. However that involves computing the product of eight $4 \times 4$ determinants, which is only possible using CAS.
The actual CAS computation to do this for $x_4$ can be summarised as follows. Compute $N$ and $D$ using \eqref{eq:key} \begin{equation*} D(x_1,x_2,x_3,y_1,y_2,y_3) x_4 - N(x_1,x_2,x_3,y_1,y_2,y_3) \space = \space \frac 1 {(x_4-x_1)(x_4-x_2)(x_4-x_3)} \NormalForm\left(\prod_{\text{both signs}} \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & x_4 & x_4^2 & \pm y_4 \\ \end{vmatrix}, \space \mathfrak{B}\right) \end{equation*} where $\mathfrak{B}$ is the Gröbner basis \begin{equation*} \mathfrak{B} \space = \space \Basis\left(\{y_i^2-R(x_i): i=1,2,3,4\}, \space \lex(y_1,y_2,y_3,y_4)\right) \end{equation*} Then to rationalise the denominator $D$, compute the rationalising factor $P$ removing unnecessary factors arising from common zeroes, \begin{equation*} P(x_1,x_2,x_3) \space = \space \frac{\NormalForm\left(D(-x_1,x_2,x_3) D(x_1,-x_2,x_3) D(x_1,x_2,-x_3), \space \mathfrak{B}'\right)} {(x_1-x_2)^2(x_1-x_3)^2(x_2-x_3)^2} \end{equation*} where $\mathfrak{B}'$ is the basis $\mathfrak{B}$ with the $x_4,y_4$ component removed. Now multiply both the numerator and denominator by $P$ and remove further common factors \begin{align*} N^*(x_1,x_2,x_3,y_1,y_2,y_3) \space &= \space \frac{\NormalForm\left(N(x_1,x_2,x_3,y_1,y_2,y_3) P(x_1,x_2,x_3), \space \mathfrak{B}'\right)} {(x_1-x_2)^2(x_1-x_3)^2(x_2-x_3)^2} \\\\ D^*(x_1,x_2,x_3) \space &= \space \frac{\NormalForm\left(D(x_1,x_2,x_3) P(x_1,x_2,x_3), \space \mathfrak{B}'\right)} {(x_1-x_2)^2(x_1-x_3)^2(x_2-x_3)^2} \\ \end{align*} We have finally \begin{equation} \label{eq:x4rat} x_4 \space = \space \frac{N^*(x_1,x_2,x_3,y_1,y_2,y_3)} {D^*(x_1,x_2,x_3)} \end{equation} There does not appear to be any simple way to express $N^*$ but we can observe that $D^*$ is essentially the level 4 addition formula evaluated at $x_4 = \infty$.
The computation for $y_4$ can be summarised as follows. Compute $\hat{N}$ and $\hat{D}$ using \begin{equation*} \hat{D}(x_1,x_2,x_3,y_1,y_2,y_3) y_4 - \hat{N}(x_1,x_2,x_3,y_1,y_2,y_3) \space = \space \frac 1 {(y_4-y_1)(y_4-y_2)(y_4-y_3)} \NormalForm\left(\prod_{i=1}^4 \begin{vmatrix} 1 & x_1 & x_1^2 & y_1 \\ 1 & x_2 & x_2^2 & y_2 \\ 1 & x_3 & x_3^2 & y_3 \\ 1 & \chi_i & \chi_i^2 & y_4 \\ \end{vmatrix}, \space \hat{\mathfrak{B}}\right) \end{equation*} where $\hat{\mathfrak{B}}$ is the Gröbner basis \begin{multline*} \hat{\mathfrak{B}} \space = \space \Basis\Big(\big\{y_1^2 - R(x_1), \space y_2^2 - R(x_2), \space y_3^2 - R(x_3), \space a(\chi_1+\chi_2+\chi_3+\chi_4) + b, \space a(\chi_1\chi_2+\chi_1\chi_3+\chi_1\chi_4+\chi_2\chi_3+\chi_2\chi_4+\chi_3\chi_4) - c, \\ a(\chi_1\chi_2\chi_3+\chi_1\chi_2\chi_4+\chi_1\chi_3\chi_4+\chi_2\chi_3\chi_4) + d, \space a\chi_1\chi_2\chi_3\chi_4 - e + y_4^2\big\}, \space \lex(x_1,x_2,x_3,\chi_1,\chi_2,\chi_3,\chi_4)\Big) \end{multline*} Then rationalise the denominator $\hat{D}$ as follows \begin{align*} \hat{N}^*(x_1,x_2,x_3,y_1,y_2,y_3) \space &= \space \frac{\NormalForm\left(\hat{N}(x_1,x_2,x_3,y_1,y_2,y_3) P^2(x_1,x_2,x_3), \space \mathfrak{B}'\right)} {(x_1-x_2)^4(x_1-x_3)^4(x_2-x_3)^4} \\\\ \hat{D}^*(x_1,x_2,x_3) \space &= \space \frac{\NormalForm\left(\hat{D}(x_1,x_2,x_3,y_1,y_2,y_3) P^2(x_1,x_2,x_3), \space \mathfrak{B}'\right)} {(x_1-x_2)^4(x_1-x_3)^4(x_2-x_3)^4} \\\\ \end{align*} As a sanity check we should get $\hat{D}^* = a^8 {D^*}^2$. We have finally \begin{equation} \label{eq:y4rat} y_4 \space = \space \frac{\hat{N}^*(x_1,x_2,x_3,y_1,y_2,y_3)} {\hat{D}^*(x_1,x_2,x_3)} \end{equation}
Lemniscatic Level 2 Addition Formula
Using CAS it is relatively easy to compute \eqref{eq:x4rat} but the resulting formula is too large to include here. Using the curve $y^2 = 1 - x^4$ we get a simpler formula which is still representative of the general formula. So if \begin{equation*} \bigint_{x_3}^{x_1} \frac {dx} {\sqrt{1-x^4}} \enspace + \enspace \bigint_{x_3}^{x_2} \frac {dx} {\sqrt{1-x^4}} \enspace = \enspace\bigint_{x_3}^{x_4} \frac {dx} {\sqrt{1-x^4}} \end{equation*} then we have \begin{equation} \label{eq:addvar4lev2lem} x_4 \space = \space \frac {2 x_1 x_2 x_3\left(x_1^2 + x_2^2 + x_3^2 - x_1^2 x_2^2 x_3^2\right) \space + \space \left(x_1^2 x_2^2 - x_1^2 x_3^2 - x_2^2 x_3^2 - 1\right)x_3 y_1 y_2 \space + \space \left(-x_1^2 x_2^2 + x_1^2 x_3^2 - x_2^2 x_3^2 - 1\right)x_2 y_1 y_3 \space + \space \left(-x_1^2 x_2^2 - x_1^2 x_3^2 + x_2^2 x_3^2 - 1\right)x_1 y_2 y_3} { \left(x_1^2 x_2^2 + x_1^2 x_3^2 + x_2^2 x_3^2 + 1\right)^2 \space - \space 4 x_2^2 x_2^2 x_3^2 \left(x_2^2 + x_2^2 + x_3^2\right)} \end{equation} and when we put $(x_3,y_3)=(0,-1)$ in \eqref{eq:addvar4lev2lem} we get Eulers[1] famous addition formula for the lemniscatic elliptic integral \begin{equation} x_4 \space = \space \frac {x_1y_2 \space + \space x_2y_1} {1 \space + \space x_1^2x_2^2} \end{equation} If we instead use the $y$ formula \eqref{eq:y4rat} for the same curve we get an addition formula for the second Schwarz-Christoffel elliptic integral \begin{equation*} \bigint_{-1}^{y_1} \frac {dy} {(1-y^2)^{\frac 3 4}} \enspace + \enspace \bigint_{-1}^{y_2} \frac {dy} {(1-y^2)^{\frac 3 4}} \enspace = \enspace \bigint_{-1}^{y_4} \frac {dy} {(1-y^2)^{\frac 3 4}} \end{equation*} then we have \begin{equation} y_4 \space = \space \frac {2 x_1 x_2 (x_1^2 \space + \space x_2^2) \space - \space (1 \space - \space x_1^2 x_2^2) y_1 y_2} {\left(1 \space + \space x_1^2 x_2^2\right)^2} \end{equation}
Addition Formulae In Three Variables
We can calculate a few more addition formula for the curve $y^2 = (1-x^2)(1-k^2x^2)$. Using \eqref{eq:x4rat} and \eqref{eq:y4rat} and putting $(x_3, y_3)$ equal to various auspicious points on the curve gives
| $(x_3,y_3)$ | $x_4$ | $y_4$ |
|---|---|---|
| $(0,-1)$ | $\displaystyle \frac {x_1 y_2 \space + \space x_2 y_1} {1 \space - \space k^2 x_1^2 x_2^2}$ | $\displaystyle \frac {(1 + k^2) x_1 x_2 (1 + k^2 x_1^2 x_2^2) - 2k^2 x_1 x_2(x_1^2 + x_2^2) - (1 + k^2 x_1^2 x_2^2) y_1 y_2} {\left(1 \space - \space k^2 x_1^2 x_2^2\right)^2} $ |
| $(1,0)$ | $\displaystyle \frac {(1-k^2) x_1x_2 \space - \space y_1y_2} {1 - k^2x_1^2 - k^2x_2^2 + k^2x_1^2x_2^2}$ | $\displaystyle \frac {(1-k^2)\big[k^2(x_1^2 - x_2^2)(x_1y_2 - x_2y_1) \space - \space (1 - k^2x_1^2x_2^2)(x_1y_2 + x_2y_1)\big]} {\left(1 - k^2x_1^2 - k^2x_2^2 + k^2x_1^2x_2^2\right)^2}$ |
| $(-\infty,\infty)$ | $\displaystyle \frac {x_1y_2 \space - \space x_2y_1} {k(x_1^2 - x_2^2)}$ | $\displaystyle \frac {\big[(1 + k^2)x_1x_2 + y_1y_2\big](x_1^2 + x_2^2) \space - \space 2x_1x_2(1 + k^2x_1^2x_2^2)} {k\left(x_1^2 - x_2^2\right)^2}$ |
Addition Formulae In Three Variables - Alternate Approach
To construct three variable addition formula we could start with the three dimensional basis $1,x,y \pm x^2\sqrt{a}$. That is the space of elliptic functions with a pole of order at most 2 at $\rho$ and order at most 1 at $-\rho$ (or vice-versa if we choose the other sign on the square root). In that case we will have a level one addition formula like \begin{equation} \label{eq:add3a} \begin{vmatrix} 1 & x_1 & y_1 \pm \sqrt{a} \thinspace x_1^2 \\ 1 & x_2 & y_2 \pm \sqrt{a} \thinspace x_2^2 \\ 1 & x_3 & y_3 \pm \sqrt{a} \thinspace x_3^2 \\ \end{vmatrix} \space = \space 0 \end{equation} However \eqref{eq:add3a} is easily seen to be equivalent to \eqref{eq:add4} when $(x_4,y_4) \rightarrow (\infty, \mp \infty)$. That is to say it is simply equivalent to the four variable addition formula with one of the points set to one of the two points on the curve at infinity, depending on the sign of the square root.
References
[1] Observationes de comparatione arcuum curvarum irrectificibilium Euler Archive - All Works. 252.