前回の投稿「楕円曲線をsympyのplot_implicitでグラフにする」では楕円曲線上の有理点の加算とスカラー倍をグラフにしました。楕円曲線暗号の仕組みに近づくためにもう一度スカラー倍をグラフを抜きにして考えます。

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

例題１: 楕円曲線$y^2=x^3+17$上の点$P(2,5)$の$2P, 4P, 8P$を求めよ

この結果を見ると$2P, 4P, 8P$は有理数点ではあるもののそれを表す数値の桁数が急激に大きくなっているのが分かりますす。このような場合の対処の仕方の常套手段として剰余をとることが考えられます。

import sympy as sp
# ----------- 2Pを計算する --------------
def mul2(a, b, pt):
  x1, y1 = pt
  l = sp.Rational((3*x1**2+a),2*y1)
  x3 = l**2-2*x1
  y3 = l*(x1-x3)-y1
  return (x3, y3) 

# 楕円曲線 y^2 = x^3 + ax + b 
a, b = 0, 17
x, y = 2, 5
for i in range(3):
  (x, y) = mul2(a, b, (x,y))
  print(2**(i+1), (x,y)
# 2 (-64/25, 59/125)
# 4 (38194304/87025, -236046706033/25672375)
# 8 (532027047589930897040873195264/4848863077511293855911670225, -388064005784387552318916270407513322740532287/337644656448214941842939018840311120390375)

剰余演算（有限体）に変更

この演算を素数$p=13$剰余で考えます。基本的には剰余計算なので加減乗はそのままpで割った余りを求めますが、割り算は注意が必要です。ここではpythonの$pow$関数でべき乗を$-1$とすること剰余の割り算を計算しています。これによって点の座標は$[0..12]$の整数に限られます。

import sympy as sp
# ----------- 2Pを計算する --------------
def mul2(a, b, pt, p):
  x1, y1 = pt
  l = (((3*x1**2+a)%p)*pow(2*y1, -1, p)) % p
  x3 = (l**2-2*x1)%p
y3 = (l*(x1-x3)-y1)%p
  return (x3, y3)

# 楕円曲線 y^2 = x^3 + ax + b
a, b = 0, 17
x, y = 2, 5
p = 13  # modulo p
for i in range(6):
  (x, y) = mul2(a, b, (x,y), p)
  print(2**(i+1), (x,y))
# 2 (12, 9)
# 4 (6, 5)
# 8 (10, 9)
# 16 (5, 5)
# 32 (4, 9)

剰余演算で点Pのスカラー倍n x Pを求める

今までは2倍するコードだけでしたが、ここからスカラー倍n x Pを求めるを作っていきます。ここでの要点は以下の３つです。

点PとQの足し算をするコードを追加
直線の傾きが無限大になった時の結果を無限遠点 Oとし、これを単位元とみなす。従って$O + P = P + O = P$となる
掛け算を高速化するために「繰り返し２乗法」を使う

加算と2倍を行う関数を作る

まず点p1とp2加算と2倍を行う関数ECaddを以下のように定義します。

無限遠を(oo,oo)とする。ooはsympy.ooで表す
もしどちらかが無限遠点だったら、もう一つの点を返す
もし傾きが無限大なら無限遠点を返す
もしp1とp2だったら2倍にする

import sympy as sp
# calculate p1 + p2 or 2 x p
def ECadd(p1, p2):
    (x1, y1), (x2, y2) = p1, p2
    if y1 == sp.oo: return p2  # p1 is infinity
    if y2 == sp.oo: return p1  # p2 is infinity
    if x2 == x1 and y1 != y2: # vertical
      return (sp.oo, sp.oo)    # return infinity
    if x2 == x1 and y1 == y2:  # same point
      if y1 == 0: return (sp.oo, sp.oo)    # return infinity
      l = (((3 * x1**2 + a) % p) * pow(2 * y1, -1, p)) % p
    else:
      l = ((y2 - y1) * pow(x2 - x1, -1, p)) % p

    x3 = (l** 2 - x1 - x2) % p
    y3 = (l * (x1 - x3) - y1) % p
    return (x3, y3)

繰り返し２乗法でn 倍を作る

繰り返し２乗法を使い、nが偶数なら2倍、奇数ならその前にそこまでの答えを足します。

# calculate n x P
def ECmul(n, pt):
    ret = (sp.oo, sp.oo)
    while n > 0:
      if n % 2 == 1:
        ret = ECadd(ret, pt)
      pt = ECadd(pt, pt)  # pt = 2 x pt
      n //= 2
    return ret

無限遠点（単位元）になるまでn倍する

楕円曲線$y^2=x^3+17$に対して$p = 13$で$P(2,5)$のn倍を計算していくと21倍で無限遠点（単位元）になりました。この後はこれが繰り返されます。

a, b = 0, 17
p = 13
p1 = (2, 5)
print(f"--- p = {p}---")
for d in range(1,2*p):
  p2 = ECmul(d, p1)
  print(d, p2)
  if p2 == (sp.oo, sp.oo):
    break
# --- p = 13---
# 1 (2, 5)
# 2 (12, 9)
# 3 (8, 3)
# 4 (6, 5)
# 5 (5, 8)
# 6 (7, 3)
# 7 (0, 2)
# 8 (10, 9)
# 9 (11, 10)
# 10 (4, 4)
# 11 (4, 9)
# 12 (11, 3)
# 13 (10, 4)
# 14 (0, 11)
# 15 (7, 10)
# 16 (5, 5)
# 17 (6, 8)
# 18 (8, 10)
# 19 (12, 4)
# 20 (2, 8)
# 21 (oo, oo)

楕円曲線上の整数点（有限体）の数

楕円曲線$y^2=x^3+17$に対して$p = 13$で$P(2,5)$のn倍を求めていくと20個の点を経て無限遠点になりました。素朴な疑問としてこれはすべての整数点（有限体）の点を網羅しているのかという疑問が出てきます。
ここでこの楕円曲線上のすべての$\pmod {13}$での整数点を求めてみます。

import sympy as sp
def ECpoints(a,b, p):
  ret = []
  for x in range(p):
    rs = (x**3 + a*x + b) % p
    ycand = sp.nthroot_mod(rs, 2, p, True)
    for y in ycand:
      ret += [(x,y)]
  return ret

a, b =  0, 17
p = 13
pts = ECpoints(a, b, p)
print(p, len(pts), pts)
# 13 20 [(0, 2), (0, 11), (2, 5), (2, 8), (4, 4), (4, 9), (5, 5), (5, 8), (6, 5), (6, 8), (7, 3), (7, 10), (8, 3), (8, 10), (10, 4), (10, 9), (11, 3), (11, 10), (12, 4), (12, 9)

整数点は20個となり上記の$P(2,5)$のn倍ですべてを網羅していることが分かりました。ではこの整数点の数と、$P$のn倍で網羅する点とはいつも一致するのでしょうか？

この辺を次回で検証していきたいと思います。