In this section we compute, in detail, various forms of the three variable addition formulae for the cubic elliptic curve.
The cubic elliptic curve is defined by the equation $F(x,y) = 0$ where
\begin{equation*}
F(x,y) \space = \space \sum_{i + j \le 3} a_{ij} \thinspace x^i y^j \space = \space a_{30} x^3 + a_{21} x^2 y + a_{12} x y^2 + a_{03} y^3 + a_{20} x^2 + a_{11} x y + a_{02} y^2 + a_{10} x + a_{01}
y + a_{00}
\end{equation*}
For convenience we also use the homogenous form of $F$
\begin{equation*}
F(x,y,z) \space = \space \sum_{i + j + k = 3} a_{ij} \thinspace x^i y^j z^k
\end{equation*}
This curve is particularly amenable to the calculation of a three variable addition formula
because it is uniformised by two order 3 elliptic functions, $x(t)$ and $y(t)$ which share the same 3 poles, say $\rho_1,\rho_2,\rho_3$.
The functions $1, x(t), y(t)$ are a basis for the three dimensional vector space of all elliptic functions with a pole of order at most 1 at $\rho_i$.
This implies, by the Extended Frobenius-Stickelberger formula, that they have an addition formula of the form
\begin{equation*}
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & x_3 & y_3 \\
\end{vmatrix} \space = \space 0
\end{equation*}
Our objective is to algebraically manipulate add3, subject to the conditions $F(x_i,y_i) = 0$ for $i=1,2,3$, into formulae of the form
\begin{equation*}
x_3 \space = \space X(x_1,y_1,x_2,y_2), \qquad y_3 \space = \space Y(x_1,y_1,x_2,y_2)
\end{equation*}
where $X$ and $Y$ are rational functions.
Symmetric Addition Formula - Simple Direct Method
Observe that add3 is linear in $x_3$ and $y_3$. Solving for $y_3$ gives a linear polynomial in $x_3$
\begin{equation*}
y_3 \space = \space \frac {x_3(y_1-y_2) \space + \space (x_1y_2 - x_2y_1)} {x_1-x_2}
\end{equation*}
If we substitute y3 into $F(x_3,y_3) = 0$ and clear the denominator using the homogenous form of $F$ we get a cubic polynomial in $x_3$.
But we already know two roots of this polynomial, namely $x_1$ and $x_2$ so we may write
\begin{equation*}
F(x_3(x_1 - x_2),\space x_3(y_1 - y_2) + (x_1y_2 - x_2y_1),\space x_1 - x_2) \space = \space (x_3 - x_1)(x_3 - x_2)(Dx_3 - N)
\end{equation*}
where $D$ and $N$ are rational functions in $x_1,x_2,y_1,y_2$ that we need to determine.
Now there is a simple trick we can use to determine $N$. Put $x_3 = 0$ in simple and solve for $N$ giving
\begin{equation*}
N \space = \space -\frac 1 {x_1 x_2} \thinspace F(0, \space x_1y_2 - x_2y_1, \space x_1 - x_2)
\end{equation*}
We can use a similar trick to determine $D$.
First divide simple by $x_3^3$ then let $x_3 \rightarrow \infty$ and solve for $D$ giving
\begin{equation*}
D \space = \space F(x_1 - x_2, \space y_1 - y_2, \space 0)
\end{equation*}
Computing $N / D$ we get a formula for $x_3$ and by a similar method a formula for $y_3$
\begin{equation*}
x_3 \space = \space -\frac {F(0, \space x_1y_2 - x_2y_1, \space x_1 - x_2)} {x_1x_2 \thinspace F(x_1 - x_2, \space y_1 - y_2, \space 0)}, \qquad
y_3 \space = \space -\frac {F(x_2y_1 - x_1y_2, \space 0, \space y_1 - y_2)} {y_1y_2 \thinspace F(x_1 - x_2, \space y_1 - y_2, \space 0)}
\end{equation*}
Thus we have succeeded in our stated goal of expressing $x_3$ and $y_3$ as rational functions of $x_1,y_1,x_2,y_2$.
Resultant Formula
Instead of the adhoc method above we can carry out the elimination of the $y$ variable using a resultant. Let
\begin{equation*}
G(x,y) \space = \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & x & y \\
\end{vmatrix}
\end{equation*}
We can infer an identity of the form
\begin{equation*}
\resultant_y(F(x,y), \space G(x,y)) \space = \space (x - x_1)(x - x_2)(Dx - N) \space + \space Q_1(x)F(x_1,y_1) \space + \space Q_2(x)F(x_2,y_2)
\end{equation*}
where $D$ and $N$ are unknown polynomials in $x_1,y_1,x_2,y_2$ and $Q_1(x)$ and $Q_2(x)$ are unknown cubic polynomials in $x$ whose coefficients are polynomials in $x_1,x_2$.
Using this resultant formula it is easily verified, by evaluating at $x=x_1$ and $x=x_2$, that
\begin{equation*}
F\left(x(x_1 - x_2), \space x(y_1 - y_2) + (x_1 y_2 - x_2 y_1), \space x_1 - x_2\right) \enspace = \enspace (x - x_1)(x - x_2)(Dx - N)
\space + \space (x - x_2)^3F(x_1,y_1) \space - \space (x - x_1)^3F(x_2,y_2)
\end{equation*}
and, by evaluating at $x=0$ and $x=\infty$, that
\begin{aligned}
N \space &= \space - \frac 1 {x_1x_2} \Big[ F(0, \space x_1y_2 - x_2y_1, \space x_1 - x_2) \space + \space x_2^3 F(x_1,y_1) \space - \space x_1^3F(x_2,y_2) \Big] \\\\
D \space &= \space F(x_1 - x_2, \space y_1 - y_2, \space 0) \space - \space F(x_1,y_1) \space + \space F(x_2,y_2)
\end{aligned}
and therefore we may write
\begin{aligned}
x \space &= \space -\frac {(x_1x_2)^{-1}\Big[ F(0, x_1y_2 - x_2y_1, x_1 - x_2) \space + \space x_2^3 F(x_1, y_1) \space - \space x_1^3 F(x_2, y_2) \Big]}
{F(x_1 - x_2, y_1 - y_2, 0) \space - \space F(x_1,y_1) \space + \space F(x_2,y_2)} \\\\
y \space &= \space -\frac {(y_1y_2)^{-1}\Big[ F(x_2y_1 - x_1y_2, 0, y_1 - y_2) \space + \space y_2^3 F(x_1, y_1) \space - \space y_1^3 F(x_2, y_2) \Big]}
{F(x_1 - x_2, y_1 - y_2, 0) \space - \space F(x_1,y_1) \space + \space F(x_2,y_2)}
\end{aligned}
where the numerator and denominator are polynomials of degree 4 and 3 respectively.
Note the $x$-formula in x3y3simple is also true for points which do not have to lie on the curve $F(x,y)=0$.
All that is required is that $x_1,x_2$ are the $x$-coordinates of the other two points of intersection.
When they do lie on the curve we get the $x$ formula in x3y3.
And the $y$ formula only requires that $y_1,y_2$ are the $y$-coordinates of the other two points of intersection.
Pythagorean Triples
Equations x3y3 can be written in homogenous form
\begin{equation*}
x_3 \space = \space \frac 1 {x_1x_2} F(0,x_1y_2 - x_2y_1,x_1z_2 - x_2z_1), \qquad
y_3 \space = \space \frac 1 {y_1y_2} F(x_2y_1 - x_1y_2,0,y_1z_2 - y_2z_1), \qquad
z_3 \space = \space \frac 1 {z_1z_2} F(x_2z_1 - x_1z_2,y_2z_1 - y_1z_2,0)
\end{equation*}
Although the numerators are degree 6 polynomials, they may be reduced modulo $\{F(x_1,y_1,z_1), F(x_2,y_2,z_2)\}$ so that they become divisible by the denominators
which in turn make the RHS's polynomials of degree 4, as follows
\begin{aligned}
x_3 \space &= \space \frac 1 {x_1x_2} \Big[ F(0,x_1y_2 - x_2y_1,x_1z_2 - x_2z_1) \space + \space x_2^3 F(x_1,y_1,z_1) \space - \space x_1^3 F(x_2,y_2,z_2) \Big] \\\\
y_3 \space &= \space \frac 1 {y_1y_2} \Big[ F(x_2y_1 - x_1y_2,0,y_1z_2 - y_2z_1) \space + \space y_2^3 F(x_1,y_1,z_1) \space - \space y_1^3 F(x_2,y_2,z_2) \Big] \\\\
z_3 \space &= \space \frac 1 {z_1z_2} \Big[ F(x_2z_1 - x_1z_2,y_2z_1 - y_1z_2,0) \space + \space z_2^3 F(x_1,y_1,z_1) \space - \space z_1^3 F(x_2,y_2,z_2) \Big]
\end{aligned}
Thus if $F$ has integer coefficients and $(x_1,y_1,z_1)$ and $(x_2,y_2,z_2)$ are two integer points on the curve $F(x,y,z) = 0$
then a third integer point $(x_3,y_3,z_3)$ is given by pythagoras
\begin{equation*}
F(x_3,y_3,z_3) \space \equiv \space 0 \mod \space \{F(x_1,y_1,z_1), \space F(x_2,y_2,z_2)\}
\end{equation*}
For example if $Ax_1^3 \space + \space By_1^3 \space + \space Cz_1^3 \space = \space 0$ and $Ax_2^3 \space + \space By_2^3 \space + \space Cz_2^3 \space = \space 0$ then
\begin{equation*}
A\Big[By_1y_2\left(x_1y_2 - x_2y_1\right) \space + \space Cz_1z_2\left(x_1z_2 - x_2z_1\right)\Big]^3\space + \space
B\Big[Ax_1x_2\left(x_2y_1 - x_1y_2\right) \space + \space Cz_1z_2\left(y_1z_2 - y_2z_1\right)\Big]^3\space + \space
C\Big[Ax_1x_2\left(x_2z_1 - x_1z_2\right) \space + \space By_1y_2\left(y_2z_1 - y_1z_2\right)\Big]^3\space = \space 0
\end{equation*}
This is a cubic analog of the Pythagorean triple formula, if $x_1^2 \space + \space y_1^2 \space - \space z_1^2 \space = \space 0$ and $x_2^2 \space + \space y_2^2 \space - \space z_2^2
\space = \space 0$ then
\begin{equation*}
\left(x_1y_2 - x_2y_1\right)^2 \space + \space \left(x_1x_2 + y_1y_2\right)^2 \space - \space \left(z_1z_2\right)^2 \space = \space 0
\end{equation*}
Determinant Product Form
Another way to express x3y3hom is
\begin{equation*}
x_3 \space = \space \frac {a_{03}} {x_1 x_2} \space \prod\limits_{i=1}^3
{\begin{vmatrix}
x_1 & y_1 & z_1 \\
x_2 & y_2 & z_2 \\
0 & \alpha_i & 1\\
\end{vmatrix}},\qquad\qquad
y_3 \space = \space \frac {a_{00}} {y_1 y_2} \space \prod\limits_{i=1}^3
{\begin{vmatrix}
x_1 & y_1 & z_1 \\
x_2 & y_2 & z_2 \\
1 & 0 & \beta_i \\
\end{vmatrix}},\qquad\qquad
z_3 \space = \space \frac {a_{30}} {z_1 z_2} \space \prod\limits_{i=1}^3
{\begin{vmatrix}
x_1 & y_1 & z_1 \\
x_2 & y_2 & z_2 \\
\gamma_i & 1 & 0\\
\end{vmatrix}}
\end{equation*}
where $\alpha_i,\beta_i,\gamma_i$ are the roots of the cubic polynomials
\begin{equation*}
F(0,y,z) = a_{03} \prod\limits_{i=1}^3 \left(y - \alpha_i z\right) \qquad\qquad
F(x,0,z) = a_{00} \prod\limits_{i=1}^3 \left(z - \beta_i x\right) \qquad\qquad
F(x,y,0) = a_{30} \prod\limits_{i=1}^3 \left(x - \gamma_i y\right)
\end{equation*}
Observe that if we expand the determinants in x3y3z3prod by their last rows then we end up with formulae x3y3hom.
The $x$ determinant product is essentially
$\resultant_{(y/z)}\Big(F(x,y,z),
\begin{vmatrix}
x_1 & y_1 & z_1 \\
x_2 & y_2 & z_2 \\
x & y & z \\
\end{vmatrix}\Big)$
evaluated at $x=0$.
The $y$ and $z$ products similarly.
Symmetric Duplication Formula
Although all the RHS's of x3y3hom vanish when $(x_2,y_2,z_2)=(x_1,y_1,z_1)$ by evaluating them along the tangent line
\begin{equation*}
(x-x_1) F_x(x_1,y_1,z_1) \space + \space (y-y_1)F_y(x_1,y_1,z_1) \space = \space 0 \qquad\qquad z \space = \space z_1
\end{equation*}
it is straight forward to calculate the homogenous limit as $(x,y,z) \rightarrow (x_1,y_1,z_1)$, giving (writing $x$ for $x_1,x_2$ and $X$ for $x_3$)
\begin{equation*}
X \space = \space \frac 1 {x^2} F(0, F_z, -F_y), \qquad
Y \space = \space \frac 1 {y^2} F(-F_z, 0, F_x), \qquad
Z \space = \space \frac 1 {z^2} F(F_y, -F_x, 0)
\end{equation*}
where $F_x,F_y,F_z$ denote the partial derivatives evaluated at $(x,y,z)$.
Further these formulae can be converted into polynomials of degree 4 like pythagoras by reducing modulo $F(x,y,z)$.
However the computations are considerably more complicated.
\begin{equation*}
X \enspace = \enspace \frac 1 {x^2} \Big[F(0,\thinspace F_z(x,y,z),-F_y(x,y,z)) \enspace - \enspace Q(x,y,z) F(x,y,z) \Big]
\end{equation*}
To compute $Q$ we seek two polynomial functions $R(y,z)$ and $S(y,z)$ such that $F(0,F_z(x,y,z),-F_y(x,y,z)) - \left(R(y,z) + xS(y,z)\right)F(x,y,z)$ is divisible by $x^2$.
There is no obvious reason why $R$ and $S$ should be polynomials rather than rational functions.
They can be found by computing the Taylor series for $F(0,F_z,-F_y)$ and $(R+xS)F$ at $x=0$ and then equating the powers of $x^0$ and $x^1$.
The Taylor series are given by
\begin{aligned}
F(0,F_z(x,y,z),-F_y(x,y,z)) \space &= \space F(0, F_z(0,y,z), -F_y(0,y,z)) \space + \space \Big[ F_{xz}(0,y,z)F_y(0, F_z(0,y,z),-F_y(0,y,z)) - F_{xy}(0,y,z)F_z(0, F_z(0,y,z),-F_y(0,y,z)) \Big]
\thinspace x + \space \bigO(x^2) \\\\
\big[ R(y,z) + x S(y,z)\big] F(x,y,z) \space &= \space R(y,z) F(0,y,z) \space + \space \big[R(y,z)F_x(0,y,z) + S(y,z)F(0,y,z)\big] x \space + \space \bigO(x^2)
\end{aligned}
Therefore $R(y,z)$ and $S(y,z)$ are given by
\begin{aligned}
R(y,z) \space &= \space \frac {F(0, F_z(0,y,z), -F_y(0,y,z))} {F(0,y,z)} \\\\
S(y,z) \space &= \space \frac {F_{xz}(0,y,z)F_y(0, F_z(0,y,z),-F_y(0,y,z)) \space - \space F_{xy}(0,y,z)F_z(0, F_z(0,y,z),-F_y(0,y,z)) \space - \space R(y,z)F_x(0,y,z)} {F(0,y,z)}
\end{aligned}
To complete the computation of $Q$ it is simply a matter of explicitly carrying out the algebraic computations on the RHS's, using CAS, to rather miraculously obtain
\begin{aligned}
R \space &= \space \left(-27a_{00}a_{03}^2 + 9a_{01}a_{02}a_{03} - 2a_{02}^3 \right)y^3 \space + \space
\left(-27a_{00}a_{02}a_{03} + 18a_{01}^2a_{03} - 3a_{01}a_{02}^2 \right)y^2z \space + \space
\left(27a_{00}a_{01}a_{03} - 18a_{00}a_{02}^2 + 3a_{01}^2a_{02} \right)yz^2 \space + \space
\left(27a_{00}^2a_{03} - 9a_{00}a_{01}a_{02} + 2a_{01}^3 \right)z^3 \\\\
S \space &= \space \left(-27a_{00}a_{03}a_{12} + 3a_{01}a_{02}a_{12} + 9a_{01}a_{03}a_{11} - 3a_{02}^2a_{11} \right)y^2 \space + \space
\left(-18a_{00}a_{02}a_{12} + 6a_{01}^2a_{12} + 18a_{01}a_{03}a_{10} - 6a_{02}^2a_{10} \right)yz \space + \space
\left(-9a_{00}a_{02}a_{11} + 27a_{00}a_{03}a_{10} + 3a_{01}^2a_{11} - 3a_{01}a_{02}a_{10} \right)z^2 \\\\
Q \space &= \space R(y,z) \space + \space x \thinspace S(y,z)
\end{aligned}
where $Q$ is a cubic polynomial given by
\begin{aligned}
Q(x,y,z) \space = \space
\left(-27a_{00}a_{03}a_{12} + 3a_{01}a_{02}a_{12} + 9a_{01}a_{03}a_{11} - 3a_{02}^2a_{11} \right)xy^2 \space + \space
\left(-18a_{00}a_{02}a_{12} + 6a_{01}^2a_{12} + 18a_{01}a_{03}a_{10} - 6a_{02}^2a_{10} \right)xyz \space + \\
\left(-9a_{00}a_{02}a_{11} + 27a_{00}a_{03}a_{10} + 3a_{01}^2a_{11} - 3a_{01}a_{02}a_{10} \right)xz^2 \space + \space
\left(-27a_{00}a_{03}^2 + 9a_{01}a_{02}a_{03} - 2a_{02}^3 \right)y^3 \space + \space
\left(-27a_{00}a_{02}a_{03} + 18a_{01}^2a_{03} - 3a_{01}a_{02}^2 \right)y^2z \space + \\
\left(27a_{00}a_{01}a_{03} - 18a_{00}a_{02}^2 + 3a_{01}^2a_{02} \right)yz^2 \space + \space
\left(27a_{00}^2a_{03} - 9a_{00}a_{01}a_{02} + 2a_{01}^3 \right)z^3 \qquad\qquad
\end{aligned}
There are similar expressions for $Y$ and $Z$.
Alternatively the reduction can be performed directly via a CAS Gröbner basis computation yielding the same formulae
\begin{aligned}
\NormalForm(X,\left[y^2 F,\space yz F, \space z^2 F\right],\lex(y,z)) \space &= \space \left(-a_{00}a_{21}^3 + a_{01}a_{20}a_{21}^2 - a_{02}a_{20}^2a_{21} + a_{03}a_{20}^3\right) x^4 \space +
\space \ldots \\
\NormalForm(Y,\left[z^2 F,\space zx F, \space x^2 F\right],\lex(z,x)) \space &= \space \left( 9 a_{00} a_{21}^2 a_{30} - 9 a_{02} a_{10} a_{30}^2 + 3 a_{02} a_{20}^2 a_{30} - 3 a_{10} a_{11}
a_{21} a_{30} + 3 a_{10} a_{12} a_{20} a_{30} - a_{10} a_{20} a_{21}^2 + a_{11} a_{20}^2 a_{21} - a_{12} a_{20}^3 \right) x^4 \space + \space \ldots \\
\NormalForm(Z,\left[x^2 F,\space xy F, \space y^2 F\right],\lex(x,y)) \space &= \space \left( 9 a_{01} a_{12} a_{30}^2 - 3 a_{01} a_{21}^2 a_{30} - 9 a_{03} a_{20}^2 a_{30} - 3 a_{10} a_{12}
a_{21} a_{30} + a_{10} a_{21}^3 + 3 a_{11} a_{12} a_{20} a_{30} - a_{11} a_{20} a_{21}^2 + a_{12} a_{20}^2 a_{21} \right) x^4 \space + \space \ldots \\
\end{aligned}
For example, for the Hesse curve $Ax^3 \space + \space By^3 \space + \space Cz^3 \space + \space Dxyz \space = \space 0$, equation duplicationreduced,
after removing a common factor of $27ABC + D^3$, yields
\begin{equation*}
X \space = \space x\left(By^3 - Cz^3\right),\qquad Y \space = \space y\left(Cz^3 - Ax^3\right),\qquad Z \space = \space z\left(Ax^3 - By^3\right)
\end{equation*}
and we have $A X^3 \space + \space B Y^3 \space + \space C Z^3 \space + \space D XYZ \space = \space 0$.
This fact is made obvious by the simple identity, true for all $x,y,z$
\begin{equation*}
A X^3 \space + \space B Y^3 \space + \space C Z^3 \space + \space D XYZ \space = \space
\left(A x^3 \space + \space B y^3 \space + \space C z^3 \space + \space D xyz \right) \left(Ax^3 - By^3\right) \left(By^3 - Cz^3\right) \left(Cz^3 - Ax^3\right)
\end{equation*}
This is the cubic analogue of the quadratic formula
\begin{equation*}
X \space = \space x^2 - y^2,\qquad Y \space = \space 2xy,\qquad Z \space = \space z^2
\end{equation*}
\begin{equation*}
X^2 \space + \space Y^2 - \space Z^2 \space = \space
\left(x^2 \space + \space y^2 \space - \space z^2\right) \left(x^2 \space + \space y^2 \space + \space z^2\right)
\end{equation*}
Interestingly cubesdup can also be written symmetrically as
\begin{equation*}
\left(A X^3 \space + \space B Y^3 \space + \space C Z^3 \space + \space D XYZ\right) \thinspace xyz \space = \space
\left(A x^3 \space + \space B y^3 \space + \space C z^3 \space + \space D xyz \right) \thinspace XYZ
\end{equation*}
Since any cubic may be linearly transformed into Hesse form it follows that the duplication formula for every cubic may be factored in this way.
However the coefficients of that linear transform and the corresponding factorisation will, in general, be in some algebraic extension of $a_{ij}$.
ANOTHER EXAMPLE
For the curve $x^2y + y^2z + xz^2 = 0$ we have $(X,Y,Z)=(-4x^2z^2 + 2xy^2z + y^4 - 2yz^3,\space -2x^3z - 4x^2y^2 + 2xyz^2 + z^4,\space x^4 + 2x^2yz - 2xy^3 - 4y^2z^2)$ and
\begin{aligned}
X^2Y + Y^2Z + XZ^2 \space = \space \big(x^2y + y^2z + xz^2 \big)\Big(18xyz(x^6+y^6+z^6) + 17(x^6y^3+y^6z^3+z^6x^3) - 36(x^3y^6+y^3z^6+z^3x^6) \\ - 3x^2y^2z^2(x^3+y^3+z^3) -
21xyz(x^3y^3+y^3z^3+z^3x^3) + 48x^3y^3z^3\Big)
\end{aligned}
Zero's And Infinities
The next two sections are the way I originally computed x3y3.
Exploiting the fact the resultant can be written as a product of $G$ evaluated at three roots of $F$ we have
\begin{equation*}
\resultant_y(F(x,y), \space G(x,y)) \space = \space a_{03} G(x,\psi_1)G(x,\psi_2)G(x,\psi_3)
\end{equation*}
where $F(x,\psi_i) = 0$.
We will need to identify the zero's and infinities of cubic.
For this purpose we define the following three cubics in homogenous form
\begin{aligned}
F(0,y,z) \space &= \space a_{03} y^3 + a_{02} y^2 z + a_{01} y z^2 + a_{00} z^3 \space = \space a_{03}(y - \alpha_1 z)(y - \alpha_2 z)(y - \alpha_3 z) \\
F(x,0,z) \space &= \space a_{30} x^3 + a_{20} x^2 z + a_{10} x z^2 + a_{00} z^3 \space = \space a_{30}(x - \beta_1 z)(x - \beta_2 z)(x - \beta_3 z) \\
F(x,y,0) \space &= \space a_{30} x^3 + a_{21} x^2 y + a_{12} x y^2 + a_{03} y^3 \space = \space a_{30}(x - \gamma_1 y)(x - \gamma_2 y)(x - \gamma_3 y)
\end{aligned}
Then the homogenous coordinates of the points where $x=0$ are $(0, \alpha_i, 1)$.
And the homogenous coordinates of the points on cubic where $y=0$ are $(\beta_i, 0, 1)$.
And the homogenous coordinates where $x=\infty,y=\infty$ ie. $z = 0$ are $(\gamma_i, 1, 0)$.
Computation Of The Determinant Product Form
To eliminate $y_3$ from add3 take the product of add3 over the three $y_3$ roots of $F(x_3,y_3) = 0$, say $\psi_1,\psi_2,\psi_3$.
Because this is essentially the resultant of $F$ and $G$ - it just differs by a multiplicative factor independent of $x_3$ and $\psi_1,\psi_2,\psi_3$ -
we may, as before, equate it to a cubic polynomial in $x_3$ with two known roots, that is
\begin{equation*}
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & x_3 & \psi_1 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & x_3 & \psi_2 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & x_3 & \psi_3 \\
\end{vmatrix}
\space = \space (x_1-x_3)(x_2-x_3)\big[N(x_1,x_2,y_1,y_2) \space - \space D(x_1,x_2,y_1,y_2) \thinspace x_3 \big]
\end{equation*}
when $F(x_1,y_1) = 0$ and $F(x_2,y_2) = 0$.
The general form of the LHS of this equation, in this context, appears in Abels proof of his addition theorem for algebraic integrals,
see Houzel pages 89,90 [2].
From this we can obtain an expression for $N$ by evaluating X3 at $x_3=0$ and an expression for $D$ by evaluating it at $x_3=\infty$ to give
\begin{equation*}
x_1 x_2 N \space = \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & 0 & \alpha_1 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & 0 & \alpha_2 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & 0 & \alpha_3 \\
\end{vmatrix}, \qquad
D \space = \space -
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
0 & 1 & \gamma_1^{-1} \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
0 & 1 & \gamma_2^{-1} \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
0 & 1 & \gamma_3^{-1} \\
\end{vmatrix}
\end{equation*}
Then after a similar computation for $y$ we obtain the following cubic addition formula
\begin{equation*}
x_3 \space = \space -\frac {\gamma_1\gamma_2\gamma_3} {x_1 x_2}
\frac {\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & 0 & \alpha_1 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & 0 & \alpha_2 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & 0 & \alpha_3 \\
\end{vmatrix}}
{\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
0 & \gamma_1 & 1 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
0 & \gamma_2 & 1 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
0 & \gamma_3 & 1 \\
\end{vmatrix}}\qquad\qquad
y_3 \space = \space -\frac {1} {y_1 y_2}
\frac {\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & \beta_1 & 0 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & \beta_2 & 0 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & \beta_3 & 0 \\
\end{vmatrix}}
{\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
0 & \gamma_1 & 1 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
0 & \gamma_2 & 1 \\
\end{vmatrix} \space
\begin{vmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
0 & \gamma_3 & 1 \\
\end{vmatrix}}
\end{equation*}
Observe that if we expand the determinants in addcubic by their last rows then we end up with formulae x3y3.
Inflection Points
It is a well known classic result [1]
that the 9 inflection points of the cubic curve $F=0$ occur precisely at the points where the Hessian vanishes
\begin{equation*}
H \space = \space
\begin{vmatrix}
F_{xx} & F_{xy} & F_{xz} \\
F_{yx} & F_{yy} & F_{yz} \\
F_{zx} & F_{zy} & F_{zz} \\
\end{vmatrix} \space = \space 0
\end{equation*}
The rough justification for this is that if in homogenous coordinates you expand $F$ in a Taylor series about a point where the Hessian vanishes,
this means there is a line through that point, along which the second degree terms vanish, and therefore that is an inflection point.
In non-homogenous coordinates, we require that the second degree terms in the Taylor expansion of $F$ vanish along the tangent line.
Or equivalently, take the limit of add3 as $x_1,x_2,x_3 \rightarrow x$ along a tangent line.
In both cases this leads to
\begin{equation*}
L \space = \space \tfrac 1 2 \left[ F_y^2F_{xx} \space - \space 2F_xF_yF_{xy} \space + \space F_x^2F_{yy} \right] \space = \space 0
\end{equation*}
Taking the limit as $x_1,x_2,x_3 \rightarrow x$ along the curve $F = 0$ in add3 gives a formula for inflection points $x,y$
\begin{equation*}
\begin{vmatrix}
1 & x & y \\
0 & -F_y & F_x \\
0 & F_yF_{xy} - F_xF_{yy} & -F_yF_{xx} + F_xF_{xy} \\
\end{vmatrix} \space = \space 0\qquad \implies \qquad F_y^2F_{xx} - 2F_xF_yF_{xy} + F_x^2F_{yy} = 0
\end{equation*}
Expressions hessian and inflection are not equivalent, but we do have this identity valid for all $x,y,z$
\begin{equation*}
z^2 H \space = \space 6 \left(F_{xx}F_{yy} \space - \space F_{xy}^2 \right) F \space - \space 8 L
\end{equation*}
Using these homogeneity identities
\begin{equation*}
xF_x + yF_y + zF_z = 3F,\qquad xF_{xx} + yF_{xy} + zF_{xz} = 2F_x, \qquad xF_{xy} + yF_{yy} + zF_{yz} = 2F_y, \qquad xF_{xz} + yF_{yz} + zF_{zz} = 2F_z
\end{equation*}
along with elementary row and column operations, expressed as matrix multiplications, gives
\begin{equation*}
z^2H \space = \left| \begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
x & y & z \\
\end{pmatrix}
\begin{pmatrix}
F_{xx} & F_{xy} & F_{xz} \\
F_{yx} & F_{yy} & F_{yz} \\
F_{zx} & F_{zy} & F_{zz} \\
\end{pmatrix}
\begin{pmatrix}
1 & 0 & x \\
0 & 1 & y \\
0 & 0 & z \\
\end{pmatrix} \right| \space = \space
\left| \begin{pmatrix}
F_{xx} & F_{xy} & F_{xz} \\
F_{yx} & F_{yy} & F_{yz} \\
2F_x & 2F_y & 2F_z \\
\end{pmatrix}
\begin{pmatrix}
1 & 0 & x \\
0 & 1 & y \\
0 & 0 & z \\
\end{pmatrix} \right| \space = \space
\begin{vmatrix}
F_{xx} & F_{xy} & 2F_x \\
F_{yx} & F_{yy} & 2F_y \\
2F_x & 2F_y & 6F \\
\end{vmatrix} \space = \space 6 \left(F_{xx}F_{yy} \space - \space F_{xy}^2 \right) F \space - \space 8 L
\end{equation*}
And thus, apart from points at infinity, the vanishing of $L$ on the curve $F=0$ does indeed coincide with the vanishing of $H$.
Substituting inflection points into the duplication formula duplication we can infer that
\begin{equation*}
H = 0 \enspace \text{and} \enspace F = 0 \quad \implies \quad x : y : z \space = \space X : Y : Z
\end{equation*}
because direct CAS computation using $X,Y$ as defined in duplicationreduced yields
\begin{equation*}
xY - yX \enspace \equiv \enspace \tfrac 1 8 F_z H \mod F
\end{equation*}
confirming that inference.
Now can we find a nice formula for the quotient ...
Crossing Points And Cusps
For a cubic curve $F(x,y)=0$ to have a crossing point at $P$ there must be two distinct tangents through $P$.
By looking at the Taylor expansion at $P$ we can see that for this to happen $F$ and both it's partial derivatives $F_x,F_y$ must vanish at $P$
for example $F = y^2 - x^3 - x^2$ at the origin.
And that means by definition that $P$ is a singular point and therefore the genus is 0 not 1.
A variation of this is for the two distinct tangents to coincide creating a cusp for example $F = y^2 - x^3$ at the origin.
Another way to get a crossing point is for $F$ to factor into distinct components for example $F = xy(x-y)$.
Then we have two or three distinct Taylor expansions at their intersection.
In both these cases a tangent line at the crossing point is a "triple" intersection point and so we expect the duplication formula at that point to collapse in the same way it collapses at an
inflection point.
Taylor Expansion Of Resultant
Let $G$ be a straight line passing through an arbitrary point $(x_1,y_1)$ with arbitrary slope given by
\begin{equation*}
G(x,y) \space = \space
\begin{vmatrix}
1 & x_1 & y_1 \\
0 & \alpha & \beta \\
1 & x & y \\
\end{vmatrix} \space = \space \alpha (y - y_1) \space - \space \beta (x - x_1)
\end{equation*}
then the Taylor series at $x = x_1$ of the resultant $\mathfrak{R}(x) = \resultant_y\left(F(x,y),G(x,y)\right)$ is
\begin{equation*}
\mathfrak{R}(x) \space = \space \alpha^3 F \space + \space \alpha^2 \left( \alpha F_x + \beta F_y \right) (x - x_1) \space + \space
\tfrac 1 2 \alpha \left(\alpha^2 F_{xx} + 2 \alpha \beta F_{xy} + \beta^2 F_{yy}\right) (x - x_1)^2 \space + \space F(\alpha, \beta, 0) (x - x_1)^3
\end{equation*}
Using this resultant formula we have
\begin{equation*}
\mathfrak{R}(x) \space = \space \alpha^3 F\left(x, \space y_1 \space + \space ({\beta}/{\alpha}) (x - x_1) \right)
\end{equation*}
By differentiating we have
\begin{aligned}
\mathfrak{R'}(x) \space &= \space \alpha^3 F_x\left(x, \space y_1 \space + \space ({\beta}/{\alpha}) (x - x_1) \right) \space + \space \alpha^2 \beta F_y(\ldots) \\
\mathfrak{R''}(x) \space &= \space \alpha^3 F_{xx}(\ldots) \space + \space 2 \alpha^2 \beta F_{xy}(\ldots) \space + \space \alpha \beta^2 F_{yy}(\ldots) \\
\mathfrak{R'''}(x) \space &= \space \alpha^3 F_{xxx}(\ldots) \space + \space 3 \alpha^2 \beta F_{xxy}(\ldots) \space + \space 3 \alpha \beta^2 F_{xyy}(\ldots) \space + \space \beta^3
F_{yyy}(\ldots) \\
\end{aligned}
and therefore
\begin{aligned}
\mathfrak{R}(x_1) \space &= \space \alpha^3 F(x_1,y_1) \\
\mathfrak{R'}(x_1) \space &= \space \alpha^3 F_x(x_1,y_1) \space + \space \alpha^2 \beta F_y(x_1,y_1) \\
\mathfrak{R''}(x_1) \space &= \space \alpha^3 F_{xx}(x_1,y_1) \space + \space 2 \alpha^2 \beta F_{xy}(x_1,y_1) \space + \space \alpha \beta^2 F_{yy}(x_1,y_1) \\
\mathfrak{R'''}(x_1) \space &= \space \alpha^3 F_{xxx}(x_1,y_1) \space + \space 3 \alpha^2 \beta F_{xxy}(x_1,y_1) \space + \space 3 \alpha \beta^2 F_{xyy}(x_1,y_1) \space + \space \beta^3
F_{yyy}(x_1,y_1)
\space = \space 6F(\alpha,\beta,0) \\
\end{aligned}
where for clarity, the $(x_1,y_1)$ arguments have been omitted from $F,F_x,F_y,F_{xx},F_{xy},F_{yy}$.
Duplication Formula Version 2
Let $(x_1,y_1)$ be a point on the curve $F = 0$ and let $G = 0$ be the equation of the tangent line at that point, then provided at least one of partial derivatives $F_x,F_y$ is non-zero
\begin{equation*}
G(x,y) \space = \space
\begin{vmatrix}
1 & x_1 & y_1 \\
0 & F_y & -F_x \\
1 & x & y \\
\end{vmatrix} \space = \space (x - x_1) F_x \space + \space (y - y_1) F_y
\end{equation*}
So using taylor putting $\alpha = F_y$ and $\beta = -F_x$ gives this formula for the resultant
\begin{equation*}
\mathfrak{R}(x) \space = \space F_y L (x - x_1)^2 \space + \space F(F_y,-F_x,0) (x - x_1)^3
\end{equation*}
where $L = \tfrac 1 2 \left(F_y^2F_{xx} - 2F_xF_yF_{xy} + F_x^2F_{yy} \right)$.
Solving $\mathfrak{R}(x)=0$ for $x$ and $G(x,y)=0$ for $y$ gives this version of the duplication formula
\begin{equation*}
x \space = \space x_1 \space - \space \frac {F_y L} {F(F_y,-F_x,0)} \qquad\qquad y \space = \space y_1 \space + \space \frac {F_x L} {F(F_y,-F_x,0)}
\end{equation*}
When $L = 0$ we have $(x,y)=(x_1,y_1)$.
This is exactly what we expect because the vanishing $L$ is the general condition for $(x_1,y_1)$ to be an inflection point of $F$,
in other words a triple point of intersection of $F$ and $G$.
Returning to the special case where both partial derivatives $F_x,F_y$ vanish we have a crossing point and there are two tangent lines at $(x_1,y_1)$ given by the two roots of
\begin{equation*}
\alpha^2 F_{xx} + 2 \alpha \beta F_{xy} \space + \space \beta^2 F_{yy} \space = \space 0
\end{equation*}
The formula for the resultant taylor reduces to
\begin{equation*}
\mathfrak{R}(x) \space = \space F(\alpha,\beta,0) (x - x_1)^3
\end{equation*}
and again, as expected, with either tangent line, we get a triple point of intersection.
Duplication Formula Version 2 In Homogenous Coordinates
Converting dup2 to homogenous coordinates equivalent to duplication gives
\begin{equation*}
X \space = \space \frac 1 {z^3} \left[x \cdot F(F_y,-F_x,0) \space - \space F_y L \right] \qquad\qquad
Y \space = \space \frac 1 {z^3} \left[y \cdot F(F_y,-F_x,0) \space + \space F_x L \right] \qquad\qquad
Z \space = \space \frac 1 {z^2} F(F_y,-F_x,0)
\end{equation*}
which implies
\begin{equation*}
xY \space - \space yX \space = \space \frac {xF_x + yF_y} {z^3} L \space \equiv \space \tfrac 1 8 F_z H \mod F
\end{equation*}
which confirms CAS computation casdup
\begin{aligned}
X \space &= \space 8 x F(F_y,-F_x,0) \space + \space z^2 F_y H \qquad\qquad\qquad Y \space = \space 8 y F(F_y,-F_x,0) \space - \space z^2 F_x H \qquad\qquad\qquad Z \space = \space 8 z
F(F_y,-F_x,0) \\\\
X \space &= \space 8 x F(-F_z,0,F_x) \space - \space y^2 F_z H \qquad\qquad\qquad Y \space = \space 8 y F(-F_z,0,-F_x) \qquad\qquad\qquad Z \space = \space 8 z F(-F_z,0,F_x) \space + \space y^2
F_x H \\\\
X \space &= \space 8 x F(0,F_z,-F_y) \qquad\qquad\qquad Y \space = \space 8 y F(0,F_z,-F_y) \space + x^2 F_z H\qquad\qquad\qquad Z \space = \space 8 z F(0,F_z,-F_y) \space - \space x^2 F_y H \\\\
\end{aligned}
Fermat Curve
The Fermat curve is defined by
\begin{equation*}
F \enspace = \enspace Ax^3 \space + \space By^3 \space + \space Cz^3
\end{equation*}
The inflection points are given by
\begin{equation*}
H \space = \space 216ABCxyz \qquad\qquad
L \space = \space 27 ABxy(A x^3 \space + \space B y^3) \qquad\qquad
z^2H \space + \space 8L \space = \space 216ABxyF
\end{equation*}
The duplication formula components, after dividing by $27ABC$, are
\begin{equation*}
X \space = \space x\left(By^3 - Cz^3\right) \qquad\qquad Y \space = \space y\left(Cz^3 - Ax^3\right) \qquad\qquad Z \space = \space z\left(Ax^3 - By^3\right)
\end{equation*}
Hesse Curve
The Hesse curve is defined by
\begin{equation*}
F \enspace = \enspace Ax^3 \space + \space By^3 \space + \space Cz^3 \space + \space Dxyz
\end{equation*}
The inflection points are given by
\begin{aligned}
H\space &= \space (216ABC + 2D^3)xyz \space - \space 6D^2 (Ax^3 + By^3 + Cz^3) \\\\
L \space &= \space 27ABxy( Ax^3 + By^3 + Dxyz ) \space - \space D^3xyz^3 \\\\
z^2 H \space + \space 8L\space &= \space \space (216ABxy - 6D^2z^2)F
\end{aligned}
The duplication formula components, after dividing by $27ABC + D^3$, are the same as the Fermat curve
\begin{equation*}
X \space = \space x\left(By^3 - Cz^3\right) \qquad\qquad Y \space = \space y\left(Cz^3 - Ax^3\right) \qquad\qquad Z \space = \space z\left(Ax^3 - By^3\right)
\end{equation*}
The addition formula components are
\begin{aligned}
X \space &= \space 3B\left(x_2y_1 - x_1y_2\right)y_1y_2 \space + \space 3C\left(x_2z_1 - x_1z_2\right)z_1z_2 \space + \space D\left(x_2^2y_1z_1 - x_1^2y_2z_2\right) \\\\
Y \space &= \space 3A\left(x_1y_2 - x_2y_1\right)x_1x_2 \space + \space 3C\left(y_2z_1 - y_1z_2\right)z_1z_2 \space + \space D\left(x_1y_2^2z_1 - x_2y_1^2z_2\right) \\\\
Z \space &= \space 3A\left(x_1z_2 - x_2z_1\right)x_1x_2 \space + \space 3B\left(y_1z_2 - y_2z_1\right)y_1y_2 \space + \space D\left(x_1y_1z_2^2 - x_2y_2z_1^2\right)
\end{aligned}
One Two Curve
The "One Two" curve is defined by
\begin{equation*}
F \enspace = \enspace Ax^2y \space + \space By^2z \space + \space Cxz^2
\end{equation*}
The inflection points are given by
\begin{aligned}
H\space &= \space 24ABCxyz - 8(A^2Cx^3 + AB^2y^3 + BC^2z^3) \\\\
L \space &= \space BC^2z^5 - 3A^3x^4y + 4AB^2y^3z^2 - 2A^2Cx^3z^2 \\\\
z^2 H \space + \space 8L\space &= \space 24A(Byz - Ax^2)F
\end{aligned}
The duplication formula components are
\begin{aligned}
X \space &= \space B(- 4AC^2x^2z^2 + 2ABCxy^2z + AB^2y^4 - 2BC^2yz^3) \\\\
Y \space &= \space C(- 4A^2Bx^2y^2 + 2ABCxyz^2 + BC^2z^4 - 2A^2Cx^3z) \\\\
Z \space &= \space A(- 4B^2Cy^2z^2 + 2ABCx^2yz + A^2Cx^4 - 2AB^2xy^3)
\end{aligned}
The addition formula components are
\begin{aligned}
X \space &= \space Ax_1x_2\left(x_2y_1 - x_1y_2\right) \space + \space 2By_1y_2\left(x_2z_1 - x_1z_2\right) \space + \space B\left(x_2y_1^2z_2 - x_1y_2^2z_1\right) \space + \space
C\left(x_2^2z_1^2 - x_1^2z_2^2\right) \\\\
\end{aligned}
Invariance Under Linear Transformations
The inhomogenous cubic curve is mapped to another inhomogenous cubic curve by a pair of linear fractional transformations in two variables (with the same denominators).
Because such transforms map straight lines to straight lines the construction of the various addition formulae is invariant.
Equivalently the homogenous cubic curve is mapped to another such curve under a three variable linear transformation, and are also invariant (commute).
For example if $\operatorname{dup}(F)(x,y,z)$ denotes the duplication operation on a point $(x,y,x)$ on a curve $F$, and $M$ is a linear transformation then we must have
\begin{equation*}
\operatorname{dup}(F \circ M)(x,y,z) \space \times \space \big( M \operatorname{dup}(F)\big) (x,y,z) \space \equiv \space 0 \mod F(x,y,z)
\end{equation*}
For the Hesse curve $F = x^3+y^3+z^3+3xyz$ we have
\begin{equation*}
\operatorname{dup}(F)(x,y,z) \space = \space \Big(x(y^3-z^3), \thinspace y(z^3-x^3), \thinspace z(x^3-y^3)\Big)
\end{equation*}
Let $M(x,y,z) \space = \space \left(x, \thinspace y, \thinspace z - x \right)$ then
\begin{equation*}
F \circ M \space = \space - 3x^2y \space + 3x^2z \space + 3xyz \space - \space 3xz^2 \space + \space y^3 \space + \space z^3
\end{equation*}
References