2017-10-22

ビットコインで使われている楕円暗号 secp256k1 をpythonで実装してopensslで確認してみる

secp256k1

ビットコインで使われている楕円暗号 secp256k1 をpythonで実装してみます。
なお、動作確認にはopensslを用います。

こちら

で示したように、計算効率を考えなければpythonで実装するのは割と容易です。

#!/usr/bin/env python
#
# secp256k1
# http://www.secg.org/SEC2-Ver-1.0.pdf
# 

# q is prime
q = 2**256 - 2**32 - 977
# l is prime
l = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141

def expmod(b,e,m):
    if e == 0: return 1
    t = expmod(b,e/2,m)**2 % m
    if e & 1: t = (t*b) % m
    return t

def inv(x):
    return expmod(x,q-2,q)

def double_pt(P):
    x = P[0]
    y = P[1]
    if y == 0: return [0, 0]
    nu = 3*expmod(x,2,q)*inv(2*y)
    x3 = expmod(nu,2,q)-2*x
    y3 = nu*(x-x3)-y
    return [x3 % q, y3 % q]

def add_pt(P, Q):
    x1 = P[0]
    y1 = P[1]
    x2 = Q[0]
    y2 = Q[1]
    if x1 == 0 and y1 == 0: return Q
    if x2 == 0 and y2 == 0: return P
    if x1 == x2:
        if (y1 + y2) % q == 0:
            return [0, 0]
        else:
            return double_pt(P)

    lm = (y1-y2)*inv(x1-x2)
    x3 = expmod(lm,2,q)-(x1+x2)
    y3 = lm*(x1-x3)-y1
    return [x3 % q, y3 % q]

def scalarmult(P, e):
    if e == 0: return [0, 0]
    Q = scalarmult(P, e/2)
    Q = add_pt(Q, Q)
    if e & 1: Q = add_pt(Q, P)
    return Q

def isoncurve(P):
    x = P[0]
    y = P[1]
    return (y**2 - x**3 - 7) % q == 0

Bx = 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798
By = 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8
B = [Bx, By]
B2 = double_pt(B)
print "q  = %x" % q
print "Bx = %x" % Bx
print "By = %x" % By
print "l  = %x" % l
if isoncurve(B):
    print "B is on curve"
else:
    assert False, "B is not on curve!"

T = scalarmult(B, l)
print "T  = (%x, %x)" % (T[0], T[1])

privkey=0x00948dda57c9964c62703b1d54f40008e351da1cc0e0a562eac4c3f7dd369c5feb
pubkey=scalarmult(B, privkey)
print "calc_pubkey = (%x, %x)" % (pubkey[0], pubkey[1])

secp256k1では、ランダムな32バイトが秘密鍵となり、秘密鍵をLittleEndianで正の整数として、ベースポイントを秘密鍵だけスカラー倍した点が公開鍵になります。
公開鍵は、x座標が32バイト、y座標が32バイトとなるので、64バイトの長さになりますが、これを、uncompress形式と呼びます。
楕円曲線が決まっているので、x座標が決まれば、y座標は2種類の値しか取らないので、x座標と２種類のうちどちらかかの情報だけでも表現できることが分かります。
これをcompress形式と呼びます。法素数は奇数と決まっているので、y座標のどちらかが偶数ならもう片方は奇数です。通常は偶数の方を取ることが多いようです。
compress形式でy座標が偶数の場合は04ではなく02となります。

opensslで鍵ペアを生成し、表示してみます。

$ openssl ecparam -name secp256k1 -genkey -out ec-priv.pem
$ openssl ec -in ec-priv.pem -text -noout
Private-Key: (256 bit)
priv:
    00:94:8d:da:57:c9:96:4c:62:70:3b:1d:54:f4:00:
    08:e3:51:da:1c:c0:e0:a5:62:ea:c4:c3:f7:dd:36:
    9c:5f:eb
pub: 
    04:39:52:76:4a:8d:90:53:26:38:53:2f:cb:7b:a0:
    b6:15:18:1e:f2:d2:2b:7d:64:a6:d3:5e:66:59:00:
    e2:42:ad:66:1b:4d:da:0a:d3:ac:24:80:ff:0f:b2:
    35:c6:22:02:43:76:4a:42:76:16:2d:36:26:b6:64:
    f4:78:c7:58:9e
ASN1 OID: secp256k1

priv:部分は32バイトあることが分かります。
pub:の先頭の04はuncompress形式であることを示しています。
続く32バイトがx座標、その後に続く32バイトがy座標を表しています。
pythonのコードではこのprivの部分をprivkeyという変数にセットしています。

privkey=0x00948dda57c9964c62703b1d54f40008e351da1cc0e0a562eac4c3f7dd369c5feb

スクリプトを実行するとこのようになります。

$ ./secp256k1.py 
q  = fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f
Bx = 79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798
By = 483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8
l  = fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141
B is on curve
T  = (0, 0)
calc_pubkey = (3952764a8d90532638532fcb7ba0b615181ef2d22b7d64a6d35e665900e242ad, 661b4dda0ad3ac2480ff0fb235c6220243764a4276162d3626b664f478c7589e)

最後のx座標とy座標が opensslコマンドで生成した鍵ペアの公開鍵の値に一致していることが分かります。

参考までにopensslのElliptic Curve実装のリンクを上げておきます。

https://github.com/openssl/openssl/tree/master/crypto/ec

openssl では多くの楕円曲線暗号に対応しています。

$ openssl ecparam -list_curves
  secp112r1 : SECG/WTLS curve over a 112 bit prime field
  secp112r2 : SECG curve over a 112 bit prime field
  secp128r1 : SECG curve over a 128 bit prime field
  secp128r2 : SECG curve over a 128 bit prime field
  secp160k1 : SECG curve over a 160 bit prime field
  secp160r1 : SECG curve over a 160 bit prime field
  secp160r2 : SECG/WTLS curve over a 160 bit prime field
  secp192k1 : SECG curve over a 192 bit prime field
  secp224k1 : SECG curve over a 224 bit prime field
  secp224r1 : NIST/SECG curve over a 224 bit prime field
  secp256k1 : SECG curve over a 256 bit prime field
  secp384r1 : NIST/SECG curve over a 384 bit prime field
  secp521r1 : NIST/SECG curve over a 521 bit prime field
  prime192v1: NIST/X9.62/SECG curve over a 192 bit prime field
  prime192v2: X9.62 curve over a 192 bit prime field
  prime192v3: X9.62 curve over a 192 bit prime field
  prime239v1: X9.62 curve over a 239 bit prime field
  prime239v2: X9.62 curve over a 239 bit prime field
  prime239v3: X9.62 curve over a 239 bit prime field
  prime256v1: X9.62/SECG curve over a 256 bit prime field
  sect113r1 : SECG curve over a 113 bit binary field
  sect113r2 : SECG curve over a 113 bit binary field
  sect131r1 : SECG/WTLS curve over a 131 bit binary field
  sect131r2 : SECG curve over a 131 bit binary field
  sect163k1 : NIST/SECG/WTLS curve over a 163 bit binary field
  sect163r1 : SECG curve over a 163 bit binary field
  sect163r2 : NIST/SECG curve over a 163 bit binary field
  sect193r1 : SECG curve over a 193 bit binary field
  sect193r2 : SECG curve over a 193 bit binary field
  sect233k1 : NIST/SECG/WTLS curve over a 233 bit binary field
  sect233r1 : NIST/SECG/WTLS curve over a 233 bit binary field
  sect239k1 : SECG curve over a 239 bit binary field
  sect283k1 : NIST/SECG curve over a 283 bit binary field
  sect283r1 : NIST/SECG curve over a 283 bit binary field
  sect409k1 : NIST/SECG curve over a 409 bit binary field
  sect409r1 : NIST/SECG curve over a 409 bit binary field
  sect571k1 : NIST/SECG curve over a 571 bit binary field
  sect571r1 : NIST/SECG curve over a 571 bit binary field
  c2pnb163v1: X9.62 curve over a 163 bit binary field
  c2pnb163v2: X9.62 curve over a 163 bit binary field
  c2pnb163v3: X9.62 curve over a 163 bit binary field
  c2pnb176v1: X9.62 curve over a 176 bit binary field
  c2tnb191v1: X9.62 curve over a 191 bit binary field
  c2tnb191v2: X9.62 curve over a 191 bit binary field
  c2tnb191v3: X9.62 curve over a 191 bit binary field
  c2pnb208w1: X9.62 curve over a 208 bit binary field
  c2tnb239v1: X9.62 curve over a 239 bit binary field
  c2tnb239v2: X9.62 curve over a 239 bit binary field
  c2tnb239v3: X9.62 curve over a 239 bit binary field
  c2pnb272w1: X9.62 curve over a 272 bit binary field
  c2pnb304w1: X9.62 curve over a 304 bit binary field
  c2tnb359v1: X9.62 curve over a 359 bit binary field
  c2pnb368w1: X9.62 curve over a 368 bit binary field
  c2tnb431r1: X9.62 curve over a 431 bit binary field
  wap-wsg-idm-ecid-wtls1: WTLS curve over a 113 bit binary field
  wap-wsg-idm-ecid-wtls3: NIST/SECG/WTLS curve over a 163 bit binary field
  wap-wsg-idm-ecid-wtls4: SECG curve over a 113 bit binary field
  wap-wsg-idm-ecid-wtls5: X9.62 curve over a 163 bit binary field
  wap-wsg-idm-ecid-wtls6: SECG/WTLS curve over a 112 bit prime field
  wap-wsg-idm-ecid-wtls7: SECG/WTLS curve over a 160 bit prime field
  wap-wsg-idm-ecid-wtls8: WTLS curve over a 112 bit prime field
  wap-wsg-idm-ecid-wtls9: WTLS curve over a 160 bit prime field
  wap-wsg-idm-ecid-wtls10: NIST/SECG/WTLS curve over a 233 bit binary field
  wap-wsg-idm-ecid-wtls11: NIST/SECG/WTLS curve over a 233 bit binary field
  wap-wsg-idm-ecid-wtls12: WTLS curvs over a 224 bit prime field
  Oakley-EC2N-3: 
    IPSec/IKE/Oakley curve #3 over a 155 bit binary field.
    Not suitable for ECDSA.
    Questionable extension field!
  Oakley-EC2N-4: 
    IPSec/IKE/Oakley curve #4 over a 185 bit binary field.
    Not suitable for ECDSA.
    Questionable extension field!

ちなみにed25519の実装はこちらにありました。
https://github.com/openssl/openssl/blob/master/crypto/ec/curve25519.c

2017-10-19

楕円関数の加法公式

暗号

楕円関数 ${y}^{2} = {x}^{3} + a {x} + b$ の形式の加法公式をメモしておきます。

アフィン座標での加算公式

1) PとQが同一点の場合。
$x_3 = {\nu}^{2} - a - x_1 - x_2$
$y_3 = {\nu} ( x_1 - x_3 ) - y_1$
${\nu} = \frac {3 {x_1}^{2} + 2 a x_1 + b} {2 y_1}$

2) $x_1 \ne x_2$ の場合
$x_3 = {\lambda}^{2} - a - x_1 - x_2$
$y_3 = {\lambda} (x_1 - x_3) - y_1$
${\lambda} = \frac {y_1 - y_2} {x_1 - x_2}$

これをアフィン座標(affine coordinate)と呼びます。

素体 $F_p$ 上の割り算は掛け算に比べて10倍〜50倍くらいの計算量なので、
割り算の回数を減らすように別の座標系を取ることが多いです。

点(x,y)の代わりに３つの値(X,Y,Z)を使い、2点(X', Y', Z'),(X, Y, Z)の間に、
X'=cX, Y'=cY, Z'=cZの関係がある場合、この2点を同一視する座標系を射影座標(projective coordinates)と呼びます。
(X, Y, Z)と(X/Z, Y/Z, 1)は同じ点ということになりますが、これを(x,y)と同一と考えます。

同じようにして、 $X'= {c}^{2} X, Y'= {c}^{3} Y, Z'= {c} Z$ の関係を使って同一視した座標系をJacobian射影座標(Jacobian projective coordinates)と呼びます。

加法公式を射影座標に変換してみると、次のようになることが簡単な計算でわかります。
1) PとQが同一点の場合
$x_3 = {\nu}^{2} - a - x_1 - x_2$
$y_3 = {\nu} ( x_1 - x_3 ) - y_1$
${\nu} = \frac {3 {x_1}^{2} + 2 a x_1 + b} {2 y_1}$

2) $x_1 \ne x_2$ の場合
$x_3 = {\lambda}^{2} - a - x_1 - x_2$
$y_3 = {\lambda} (x_1 - x_3) - y_1$
${\lambda} = \frac {y_1 - y_2} {x_1 - x_2}$

Projective座標系での加算公式

式は単純に $x_1 = X_1 / Z_1, y_1 = Y_1 / Z_1$ を代入すれば、得られます。
割り算を全てZ座標に集めることによって、割り算の数をProjective座標とアフィン座標の変換時のみにすることができます。

1) PとQが同一点の場合
$w = z {Z_1}^{2} + 3 {X_1}^{2}$
$s = Y_1 Z_1$
$B = X_1 Y_1 a$
$h = {w}^{2} - 8B$

$X_3 \to 2 h s$
$Y_3 \to w ( 4B - h ) - 8 {Y_1}^{2} {s}^{2}$
$Z_3 \to 8 {s}^{3}$

2) $x_1 \ne x_2$ の場合
$u = Y_2 Z_1 - Y_1 Z_2$
$v = X_2 Z_1 - X_1 Z_2$
$A = {u}^{2} Z_1 Z_2 - {v}^{3} - 2 {v}^2 X_1 Z_2$

$X_3 \to v A$
$Y_3 \to u ( {v}^{2} X_1 Z_2 - A) - {v}^{3} Y_1 Z_2$
$Z_3 \to {v}^{3} Z_1 Z_2$

参考:

https://written.4403.biz/source/ecc_rev30r.pdf

https://link.springer.com/content/pdf/10.1007%2F3-540-49649-1_6.pdf

2017-10-18

secp256k1仕様

secp256k1

secp256k1に関するメモです。後ろから二番目のkは数学者Neal Koblitzのkらしいです。

曲線

${y}^2 = {x}^3 + 7$

法素数

$p = {2}^{256} - {2}^{32} - 977$
= 0xFFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE FFFFFC2F
= 115792089237316195423570985008687907853269984665640564039457584007908834671663
256ビット、10進で78桁

ベースポイント:

Bx = 0x79BE667E F9DCBBAC 55A06295 CE870B07 029BFCDB 2DCE28D9 59F2815B 16F81798
= 55066263022277343669578718895168534326250603453777594175500187360389116729240
256ビット、10進で77桁

By = 0x483ADA77 26A3C465 5DA4FBFC 0E1108A8 FD17B448 A6855419 9C47D08F FB10D4B8
= 32670510020758816978083085130507043184471273380659243275938904335757337482424
256ビット、10進で77桁

ベースポイントの位数(=素数)

L = 0xFFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE BAAEDCE6 AF48A03B BFD25E8C D0364141
= 115792089237316195423570985008687907852837564279074904382605163141518161494337
256ビット、10進で78桁

ゼロ点は無限遠点(O)となる。
python実装上は、無限遠点はx=0, y=0で表現する。

BをL倍すると最初に無限遠点になる。
$L * G = O$
コファクターは1なので、Lがこの曲線上の有理点の数(無限遠点を含む)であり、巡回群の位数である。