はじめに（Introduction）

Ed25519とはエドワーズ曲線デジタル署名アルゴリズム（Edwards-curve Digital Signature Algorithm）を実装したものです。
楕円曲線暗号の一つで、他のものより高速でコンパクトである事が特徴です。

ここでは、仕様をみながら基本的な実装していきます。

Ed25519

EdDSAパラメータ(EdDSA parameters)

EdDSAには7つのパラメーターがあります。
Ed25519のパラメータも共に記載します。

EdDSA	Ed25519
整数 $b \geq 10$	$b = 256$
$2b$ビット出力する暗号化ハッシュ関数 $\text{H}$	$\text{H}$ は SHA-512
$1 \pmod{4}$ となる素数 $q$	$q$ は素数 $2^{255} − 19$
有限体 $F_q$ 要素の $(b-1)$ ビットエンコーディング	エンコーディングはリトルエンディアン
$F_q$ の非平方要素 $d$	$d = −\frac{121665}{121666} \in F_q$
※を満たす $2^{b-4}$ と $2^{b-3}$ の間の素数 $\ell$	$\ell$は素数 $2^{252} + 27742317777372353535851937790883648493$
集合の要素 $B\ne (0,1)$	$x$ が正の一意の点 $(x, \frac{4}{5}) \in E$

以下の楕円曲線を使います。

E = \{ (x, y) \in F_q \times F_q : −x^2 + y^2 = 1 + dx^2y^2 \}

楕円曲線上の点に対して、以下の加算則を定めると、中立元（neutral element）$0=(0,1)$ で群をなします。

(x_1, y_1) + (x_2, y_2) = \left( \frac{x_1y_2 + x_2y_1}{1 + dx_1x_2y_1y_2}\ ,\ \frac{y_1y_2 + x_1x_2}{1 − dx_1x_2y_1y_2} \right)

パラメータの**※**は、$\ell B=0$ を満たす事です。（$nB$ は群における $B$ の $n$ 倍を意味します。）

要素 $(x,y)\in E$ を、$256$ビット（$b$ビット）でエンコードします。
$y$ は $255$ビット（$(b-1)$ビット）でエンコードされ、最後のビットは符号ビットです。
$x$ が負の場合、符号ビットは $1$ です。

import hashlib


class Ed25519P():
    p = 2**255 - 19
    L = 2**252 + 27742317777372353535851937790883648493
    d = (-121665 * pow(121666, p-2, p)) % p

    def __init__(self, x, y):
        if x is None or y is None:
            raise ValueError("x or y are None")
        self.x = x % self.p
        self.y = y % self.p
        # : −x^2 + y^2 = 1 + dx^2y^2
        if (-1 * pow(self.x, 2, self.p) + pow(self.y, 2, self.p)) % self.p != (1 + self.d * pow(self.x, 2, self.p) * pow(self.y, 2, self.p)) % self.p:
            raise ValueError("x and y are not on the curve")

    def __add__(self, other):
        x1 = self.x
        y1 = self.y
        x2 = other.x
        y2 = other.y
        #        x1y2 + x2y1
        # x3 = ---------------
        #       1 + dx1x2y1y2
        x3 = ((x1*y2 + x2*y1) * pow((1 + self.d*x1*x2*y1*y2) %
                                    self.p, self.p-2, self.p)) % self.p
        #        y1y2 + x1x2
        # y3 = ---------------
        #       1 − dx1x2y1y2
        y3 = ((y1*y2 + x1*x2) * pow((1 - self.d*x1*x2*y1*y2) %
                                    self.p, self.p-2, self.p)) % self.p
        return self.__class__(x3, y3)

    def __sub__(self, other):
        return self + other.__class__(self.p - self.x, self.y)

    def __rmul__(self, coefficient):
        coef = coefficient % self.L
        current = self.__class__(self.x, self.y)
        result = self.__class__(0, 1)
        while coef:
            if coef & 1:
                result += current
            current += current
            coef >>= 1
        return result

    def __eq__(self, other):
        return (self.x == other.x) and (self.y == other.y)

    def __str__(self):
        return '({:064x},{:064x})'.format(self.x, self.y)

    def encode(self):
        return int.to_bytes(self.y | ((self.x & 1) << 255), 32, 'little')

コードについて、いくつか解説します。
除算に関しては、乗法逆元の乗算としています。（参照：付録（Appendix）- 乗法逆元（multiplicative inverse））
減算（__sub__）の部分は、点を反転（正負を逆転）して加算しています。
点の反転は、点 $(x_1,y_1)$ の $x_1$ を反転したもの点$(-x_1,y_1)$ です。（参照：付録（Appendix）- 符号反転（negation））
スカラー倍算（__rmul__）の部分は二進展開を使います。（参照：付録（Appendix）- 二進展開（binary expansion））

$y$ と符号ビットから、以下の式で $x$ を求める事ができます。
$x=\pm\sqrt{(y^2-1)/(dy^2+1)}$

def DecodePoint(encode_point):
    if len(encode_point) != 32:
        raise ValueError("invalid length of encode point")
    y = int.from_bytes(encode_point, 'little')
    sign = y >> 255
    y &= (1 << 255) - 1
    # x^2 = (y^2−1)/(dy^2+1)
    x2 = ((y*y-1)
          * pow(Ed25519P.d*y*y+1, Ed25519P.p-2, Ed25519P.p)) % Ed25519P.p
    if x2 == 0:
        if sign:
            raise ValueError("invalid sign of encode point")
    # x = ±√(y^2−1)/(dy^2+1)
    x = pow(x2, (Ed25519P.p+3) // 8, Ed25519P.p)
    if (x*x - x2) % Ed25519P.p != 0:
        x = (x * pow(2, (Ed25519P.p-1) // 4, Ed25519P.p)) % Ed25519P.p
        if (x*x - x2) % Ed25519P.p != 0:
            raise ValueError("invalid square of encode point")
    if (x & 1) != sign:
        x = Ed25519P.p - x
    return Ed25519P(x, y)

By = (4*pow(5, Ed25519P.p-2, Ed25519P.p)) % Ed25519P.p
B = DecodePoint(By.to_bytes(32, 'little'))

コードについて、いくつか解説します。
パラメータの長さ $32$ バイト（$256$ビット）をチェックします。
リトルエンディアンで数値化した後、$y$ と符号（偶数or奇数）を取得します。
$y$ から、$x^2$ を計算し、平方根を求めます。（参照：付録（Appendix）- 平方根（square root））
符号（偶数or奇数）から、値を確定します。

EdDSAの鍵と署名検証(EdDSA keys and signatures)

秘密鍵と公開鍵（secret key and public key）

EdDSA _秘密鍵__は $32$ バイト（$256$ビット）のデータ $k$ です。
これをハッシュした値 $\text{H}(k)=(h_0,h_1,\ldots ,h{511})$ から数値を求めます。

a = 2^{254}+\sum_{3 \leq i \leq 253}2^ih_i\in\{ 2^{254},2^{254}+8,\ldots ,2^{255}-8 \}

捕捉します。
$254$ビット目を $1$ 、最初の $3$ビットを $0$ としたものをリトルエンディアンで数値にしたものです。

0	1	2	3	4	$\cdots$	253	254	255
$0$	$0$	$0$	$h_3$	$h_4$	$\cdots$	$h_{253}$	$1$	$0$

したがって、最小値は $2^{254}$ 、下位 $3$ビットは $0$ なので、$8$ 刻みの数値となり、最大値は $2^{255}-8$ となります。

$a$ から、$A=aB$ を計算します。
EdDSA __公開鍵__は $A$ です。

署名（signature）

__秘密鍵 $k$によるメッセージ $M$ の署名は、次のように定義されます。
秘密鍵 $k$ から求めたハッシュ値の後半部分とメッセージ $M$ を使って以下を求めます。
$r = \text{H}(h{256},\ldots ,h{511},M) \in \{0,1,\ldots , 2^{512}-1 \}$
ハッシュ値をリトルエンディアンで数値にして、$\{0,1,\ldots , 2^{512}-1 \}$ の整数とます。
$r$ から、$R=rB$ を計算します。
$a,A,r,R,M$ を使って以下を求めます。
$S = (r + \text{H}(R, A, M)a) \mod \ell$

秘密鍵 $k$ によるメッセージ $M$ の署名は、$64$ バイト（$512$ビット）の $(R,S)$ となります。
$S$ は、$32$ バイト（$256$ビット）リトルエンディアンエンコーディングです。

def Sign(secret, m):
    if len(secret) != 32:
        raise ValueError("invalid length of secret")
    h = hashlib.sha512(secret).digest()
    b = (h[0] & 0xf8).to_bytes(1, byteorder='little') + \
        h[1:31] + \
        ((h[31] & 0x7f) | 0x40).to_bytes(1, byteorder='little')
    s = int.from_bytes(b, 'little')
    A = s*B
    prefix = h[32:]
    r = int.from_bytes(hashlib.sha512(prefix + m).digest(), 'little')
    R = r*B
    k = int.from_bytes(hashlib.sha512(
        R.encode() + A.encode() + m).digest(), 'little')
    S = (r + k * s) % Ed25519P.L
    return R.encode() + S.to_bytes(32, 'little')

検証（verication）

__公開鍵 $A$__によるメッセージ $M$ と署名 $(R,S)$ の検証は、次のように行います。
$A$ は $A \in E$ として、$(R,S)$ は解析して $R \in E$ と $S \in \{0,1,\ldots , \ell-1\}$ にします。
検証者は $\text{H}(R,A,M)$ を計算し、以下をチェックします。

8SB \overset{?}{=}8R + 8\text{H}(R, A, M)A

def Verify(A, signature, m):
    if len(signature) != 64:
        raise ValueError("invalid length of signature")
    R = DecodePoint(signature[:32])
    S = int.from_bytes(signature[32:], 'little')
    H = int.from_bytes(hashlib.sha512(
        R.encode() + A.encode() + m).digest(), 'little')
    return 8*S*B == 8*R + 8*H*A

おわりに（Conclusion）

単純な実装はできましたが、さらに高速に演算することができます。
また、公開鍵や秘密鍵を $32$ バイト単位にする事によって、CPUの演算的にも高速に演算可能となっています。
それにより、コンパクトにもなっています。

個人的には、署名がSchnorrなので署名の加算が可能である点、署名のロジックにナンス生成も含まれている事から、かなり使いやすい実装だと感じました。

参考（References）

High-speed high-security signatures
http://ed25519.cr.yp.to/ed25519-20110926.pdf
RFC 8032 - Edwards-Curve Digital Signature Algorithm (EdDSA)
https://tools.ietf.org/html/rfc8032
Quadratic residue - Wikipedia
https://en.wikipedia.org/wiki/Quadratic_residue#Prime_or_prime_power_modulus

付録（Appendix）

乗法逆元（multiplicative inverse）

有限体における乗法の逆元、いわゆる $a$ に対して、$\frac{1}{a}$ となる値です。
フェルマーの小定理 $a^{p-1} \equiv 1 \pmod{p}$ より、

a^{p-2} \equiv a^{-1} \pmod{p}

で求める事ができます。

二進展開（binary expansion）

ポイントの加算において、$nB$ のように $n$倍する場合、単純には $n$回 $B$ を加算すればいいのですが、
$n$ が大きな数の場合（$2^{252} \leq \ell \leq 2^{253}$）、加算に時間がかかってしまいます。
$n$ を二進数に展開し、$B$ を倍々にし $1$ のところだけを加算する事で $\log_2(n)$回ですます事ができます。

例えば、$100$を考えてみます。

\begin{align*}
100_{10} &= 1100100_2 \\
&= 1000000_2 \\
&+ \ \ 100000_2 \\
&+ \ \ \ \ \ \ \ \ 100_2 \\
\end{align*}

通常は $100$ 回加算する必要がありますが、この方法を使えば
二倍の計算が $6$ 回、加算が $3$ 回の合計 $9$ 回で済ます事ができます。

符号反転（negation）

点 $(x_1,y_1)$ と点$(-x_1,y_1)$ を加算すると、点 $(0,1)$ となるか見てみます。
点 $(x_1,y_1)$ が楕円曲線上の点である式と加算の式で求められます。

\begin{align*}
−x_1^2 + y_1^2 &= 1 + dx_1^2y_1^2 \\
(x_1, y_1) + (-x_1, y_1) &= \left( \frac{x_1y_1 - x_1y_1}{1 - dx_1^2y_1^2}\ ,\ \frac{y_1^2 - x_1^2}{1 + dx_1^2y_1^2} \right) \\
&= \left( 0\ ,\ \frac{y_1^2 - x_1^2}{1 + dx_1^2y_1^2} \right) \\
&= \left( 0\ ,\ 1 \right) \\
\end{align*}

平方根（square root）

有限体 $F_p$ での $p \equiv 5 \pmod{8}$ の場合、$a$ が平方数である場合（$a=x^2$）、$x$ は以下の式で求める事ができます。

x =
\begin{cases}
\pm a^{\frac{p+3}{8}}\\
\pm a^{\frac{p+3}{8}}2^{\frac{p-1}{4}}\\
\end{cases}

数学的な証明については割愛します。（書ききれないTT)

プログラムとしてのロジックは、まず平方数 $a$ の $\frac{p+3}{8}$乗、すなわち $a^{\frac{p+3}{8}}$ を計算します。
これを二乗して、元の平方数 $a$ との差が$0$であれば、その値となります。
$0$ でない場合、算出した値に$2^{\frac{p-1}{4}}$を乗算します。
再びこれを二乗して、元の平方数 $a$ との差が$0$であれば、その値となります。
$0$ でない場合、$a$ は平方数ではありません。