前回ではsympyのplot_implicitが生成した座標を取り出してmatplotlibで描く方法を紹介しました。今回はそれを使って楕円曲線暗号の基礎の基礎となる2点の和と2倍したときの点の座標を描いてみます。

楕円曲線上の2点の和

例題１: 楕円曲線$y^2=x^3+17$上の2点$P_1(-1,4), P_2(2,5)$の和を求めてグラフに描け

楕円曲線暗号 (Wikipedia) の楕円曲線の加算公式より

\begin{align}
&楕円曲線 y^2 = x^3+ax+b 上の\\
&2点P_1=(x_1,y_1), P_2=(x_2,y_2)\\
&を通る直線が楕円曲線と交わる点をP_1*P_2= (x_3, y_3)とすると\\
&P_1+P_2 = (x_3,-y_3)となる \\ \\
&直線の式をy = λx+νとすると\\
&λ= (y_2-y_1)/(x_2-x_1) \\
&ν= y_1-λx_1=y_2-λx_2 \\
&これからx_3, y_3は以下のように求まる \\ \\
&x_3=λ^2-x_1-x_2 \\
&y_3 = λ(x_3-x_1)+y_1
\end{align}

これをそのままコードに下のがこれです。キャンバスの定義と楕円曲線の座標の取得は前回の投稿のinit_subplotとplot_fxyを使っています。

import sympy as sp
import matplotlib.pyplot as plt

def draw_line(ax, p1, a, CY): # p1上の傾きaの直線をキャンバスいっぱいに描く 
  (x, y), (cb, ct) = p1, CY
  ax.plot([x+(cb-y)/a, x+(ct-y)/a], CY,  lw=0.5, color='green')
  return

def plot_dot(x,y,dx,dy,pn):
  ax.plot(x,y, marker='o',  ms=3,color='blue')
  ax.text(x+dx,y+dy, pn+f"({x},{y})", fontsize="xx-small")
  return

x, y = sp.symbols("x y")
# キャンバスの定義
sizex = 6
CX, CY = (-sizex, sizex), (-sizex, sizex)
ax = init_subplot(plt, CX, CY)
plt.rcParams["font.size"] = 12
# 楕円曲線を描く
a, b = 0, 17
ax.fill(*plot_fxy(x**3+a*x+b-y**2, CX, CY), facecolor='blue')

x1, y1, x2, y2 = -1, 4, 2, 5
plot_dot(x1,y1,-2.5,0.5, r"$P_1$")  #  Plot P1
plot_dot(x2,y2,1,-0.5, r"$P_2$")    #  Plot P2
# ----------- P*Qを計算する --------------
l = sp.Rational(y2-y1,x2-x1)
x3 = l**2-x1-x2
y3 = l*(x3-x1)+y1

draw_line(ax, (x1,y1), l, CY)
plot_dot(x3,y3,0.5,-0.5, r"$P_1*P_2$")    #  Plot P1*P2
plot_dot(x3,-y3,0.5,-0.5, r"$P_1+P_2$")    #  Plot P1+P2
draw_line(ax, (x3,y3), sp.oo, CY)
plt.show

楕円曲線上の点Pから2Pを求める

例題１: 楕円曲線$y^2=x^3+17$上の点が$P(-1,4)$のとき$2P$をグラフ上に描け

これはPとQが重なったと考えればよいので傾き$\lambda$を微分で求めるだけが違いとなります。

\begin{align}
&楕円曲線 y^2 = x^3+ax+b 上の\\
&点P_1=(x_1,y_1)で接する直線\\
&楕円曲線と交わる点をP_1*P_2= (x_3, y_3)とする\\
&直線の式をy = λx+νとすると\\
&λ = (3x^2+a)/2y  \\
&ν= y_1-λx_1=y_2-λx_2 \\
&これからx_3, y_3は以下のように求まる\\ \\
&x_3=λ^2-2x_1\\
&y_3 = λ(x_3-x1)+ν
\end{align}

したがって上記のコードの$P_2$関連を削除し「P*Qを計算する」部分を以下に変えればよいですね。

# ----------- P*Qを計算する --------------
l = sp.Rational((3*x1**2+a),2*y1)
x3 = l**2-2*x1
y3 = l*(x3-x1)+y1

（開発環境：Google Colab）