La fonction $\wp$ de Weierstrass vérifie l'équation différentielle
$$ [\wp '(z)]^{2}=4[\wp (z)]^{3}-g_{2}\wp (z)-g_{3}$$de sorte que $(x,y)=(\wp(t),\wp'(t))$ fournit un paramétrage de la courbe elliptique d'équation $$y^2=4x^3-g_2x-g_3$$.
Par ailleurs, elle vérifie la formule d'addition
$$ \wp (z+y)={\frac {1}{4}}\left\{{\frac {\wp '(z)-\wp '(y)}{\wp (z)-\wp (y)}}\right\}^{2}-\wp (z)-\wp (y) $$et la formula de duplication
$${\displaystyle \wp (2z)={\frac {1}{4}}\left\{{\frac {\wp ''(z)}{\wp '(z)}}\right\}^{2}-2\wp (z),}$$valable si $2z$ n'est pas une période.
Ces formules proviennent des faits suivants : (1) toute droite coupe la courbe en 3 points (si on compte les points de tangence avec multiplicité 2), et (2) la somme des zéros d'une fonction elliptique quelconque dans un parallélograme fondamental est congrue à la somme des pôles modulo le réseau des périodes.
Ainsi, si $y=px+q$ est l'équation de la droite passant par $P_1$ et $P_2$ de paramètres $t_1$ et $t_2$, et si $t_3$ est le paramètre du troisième point d'intersection $P_3$, alors $t_1+t_2+t_3\equiv 0 \mod (\omega_1,\omega_2)$. En effet, la seule singularité de la fonction $\wp'(z)-p\wp(z)-q$ est un pôle triple en 0, la somme des trois zéros est donc un élément du réseau.
Donc, dans le groupe ${\mathbb C}/({\mathbb Z}\omega_1\oplus{\mathbb Z}\omega_2)$, on a $t_3 = -(t_1+t_2)$. Comme $\wp$ est paire, et sa dérivée $\wp'$ impaire, le point de paramètre $-t_3$ est $(x_3,-y_3)$, donc, le symétrique de $P_3$ par rapport à l'axe des abscisses.
C'est ce point que nous allons maintenant calculer pour la courbe $$y^2 = x^3+ax+b.$$
from sympy import *
var('x y x1 y1 x2 y2 s a b')
(x, y, x1, y1, x2, y2, s, a, b)
s = (y2-y1)/(x2-x1); s
E = y**2-x**3-a*x-b
L'équation de la droite $P_1P_2$ est $$\left|\begin{matrix}x-x_1& x_2-x_1\\ y-y_1&y_2-y1\end{matrix}\right|=0$$ ou encore $$y = s(x-x_1)+y_1.$$ On substitue cette valeur de $y$ dans l'équation de la courbe :
E.subs({y:(x-x1)*s+y1})
simplify(_*(x2-x1)**2)
expand(_)
On utilise ensuite le fait que $(x_1,y_1)$ vérifie l'équation de la courbe :
_.subs({y1**2:x1**3+a*x1+b})
expand(_)
factor(_)
On peut maintenant diviser par $x-x_1$ :
simplify(_/(x-x1))
On recommence avec $(x_2,y_2)$ :
_.subs({y2**2:x2**3+a*x2+b})
factor(_)
simplify(_/(x-x2))
Il ne reste plus qu'à extraire la valeur de $x_3$ :
solve(_,x)
[(a*x1 + a*x2 + 2*b + x1**2*x2 + x1*x2**2 - 2*y1*y2)/(x1**2 - 2*x1*x2 + x2**2)]
factor(_[0])
En remplaçant $ax_1+b$ par $y_1^2-x_1^3$ et $ax_2+b$ par $y_2^2-x_2^3$, on voit que cette expression vaut $$x_3=s^2-x_1-x_2.$$ Et pour finir, $$y_3 = -s(x_3-x_1)-y_1$$ d'après l'équation de la droite $P_1P_2$.
Si $P_2=P_1$, la droite $P_1P_2$ est remplacée par la tangente, qui a pour équation $$\frac{\partial F}{\partial x}(x-x_1)+\frac{\partial F}{\partial y}(y-y_1)=0$$ avec $F(x,y)=y^2-x^3-ax-b$, soit $$y = s(x-x1)+y_1$$ avec cette fois $$s=\frac{3x_1^2+a}{2y_1}.$$