はじめに

elliptic curve/secp256k1（楕円曲線）上にある点の加法をbiginteger（任意精度整数）を使って実装しみようという試みです。
ライブラリを使った方が良いですが、原理を知る事も重要なのでやってみようと思います。
※：ライブラリの方が計算速度は速いです。

Bitcoinでは、secp256k1を採用しているので、これを使います。

楕円曲線（elliptic curve）

楕円曲線は次の方程式となります。

\begin{align*}
y^2=x^3+ax+b
\end{align*}

加法

楕円曲線上の点どうしの加法を定義します。

加算：$\ R = P + Q$
楕円曲線上の２点、$P:(x_P,y_P),Q:(x_Q,y_Q)$を通る直線が楕円曲線と交わる点を$-R:(x_R,y_R)$とし、$x$軸に対象な点を$R:(x_R,-y_R)$とします。
2倍算：$\ R = 2P$
また、点$\ P\ $と点$\ Q$が同じ点であった場合、その接線と交わる点を$-R:(x_R,y_R)$とし、$x$軸に対象な点を$R:(x_R,-y_R)$とします。

楕円曲線 - Wikipediaによると以下のように計算できるそうです。

$P:(x_P,y_P),Q:(x_Q,y_Q),R:(x_R,-y_R)$とした時、加算では以下となります。

\begin{align*}
s &= \frac{y_P-y_Q}{x_P-x_Q} \\
x_R &= s^2 - x_P - x_Q \\ 
y_R &= y_P + s(x_R - x_P) \\
\end{align*}

2倍算は以下となります。

\begin{align*}
s &= \frac{3x_P^2+a}{2y_P} \\
x_R &= s^2 - 2x_P \\ 
y_R &= y_P + s(x_R - x_P) \\
\end{align*}

結合律

加法を定義しました、これらは結合律が成り立ちます。
まだ、無限遠点$O$を定義します。

\begin{align*}
P + Q &= Q + P = R \\
R - R &= O \\
P + O &= O + P = P \\
\end{align*}

なんとなく、グラフを見ているとわかると思います。

なんちゃって計算

なんとなく高校数学レベルで加法の式が求められると思いやってみましたが、案外大変でした。
とりあえず、直感的に理解出来る程度にまとめてみました。
※：偏微分って高校数学でやったっけ？

以下の楕円曲線と直線で求められる３点を考えます。

\begin{align*}
y^2 &= x^3 + ax + b \\
y &= sx + t \\
P &:(x_P,y_P) \\
Q &:(x_Q,y_Q) \\
R &:(x_R,-y_R) \\
R &= P + Q \\
\end{align*}

傾き$s$を求めます、$x_p\ne x_Q$の場合は、以下のとおりです。

\begin{align*}
y_p &= sx_p + t \\
y_Q &= sx_Q + t \\
s &= \frac{y_p-y_Q}{x_p-x_Q} \\
\end{align*}

$x_p= x_Q$の場合は、接線となるので傾きは偏微分で求めます。

\begin{align*}
y^2 &= x^3 + ax + b \\
f(x,y) &= y^2 - x^3 - ax - b = 0 \\
f'_y &= \frac{\partial f(x,y)}{\partial y} = 2y \\
f'_x &= \frac{\partial f(x,y)}{\partial x} = - 3x^2 - a \\
s &= -\frac{f'x}{f'y} = \frac{3x^2+a}{2y} \\
&x_P = x_Q \\
s &= \frac{3x_P^2+a}{2y_P} \\
\end{align*}

楕円曲線と直線から$y$を消します。

\begin{align*}
y^2 &= x^3 + ax + b \\
y &= sx + t \\
s^2x^2+2stx+t^2 &= x^3 + ax + b \\
0 &= x^3 - s^2x^2  + (a-2st)x + b-t^2 \\
\end{align*}

この式は３点を通ることから、以下の式と同等になります。

\begin{align*}
0 &= (x-x_P)(x-x_Q)(x-x_R) \\
0 &= x^3-(x_P+x_Q+x_R)x^2+(x_Px_Q+x_Qx_R+x_Rx_P)x-x_Px_Qx_R\\
\end{align*}

$x^2$の項が同じになるので、以下のように$x_R$が求まります。

\begin{align*}
s^2 &= x_P+x_Q+x_R \\
x_R &= s^2 - x_P - x_Q \\
\end{align*}

$y_R$は、直線を使って求めます。

\begin{align*}
y_R &= sx_R + t \\
-)\quad \quad y_P &= sx_P + t \\
\hline
y_R-y_P &= s(x_R-x_P) \\
y_R &= y_P + s(x_R-x_P) \\
\end{align*}

という事で楕円曲線上の点における加法が求まります。

剰余演算（modulo）

コンピュータでは桁数を無限には使えません。
そこで、剰余演算を使います。
簡単に言うと、ある数で割った余りで計算する事です。
結果、ある一定の長さ以下で表現が可能となります。

剰余演算 - Wikipedia

フェルマーの小定理

$\ p$が素数の場合、以下の式が成り立ちます。

\begin{align*}
a^{p-1} \equiv 1 \mod{p}
\end{align*}

証明は比較的簡単です。
「$a^p \equiv a \mod{p}$」が成り立つと仮定します。
「$(a+1)^p$」を二項定理で展開します。これを素数$p$で割ると最初と最後の項以外割りきれます。
仮定より、「$a+1$」でも成り立つ事がわかります。
「$a=1$」の時、「$1^p \equiv 1 \mod{p}$」で成り立つ為、帰納法により成り立つ事が証明出来ます。

\begin{align*}
(a+1)^p &= \sum_{k=0}^{p} \binom{p}{k}a^{p-k} \\
&= a^p + 1 \mod{p} \\
&= a + 1 \mod{p} \\
\end{align*}

二項定理 - Wikipedia
フェルマーの小定理 - Wikipedia

逆数（reciprocal）

$\ a^{-1}$は$s$を計算する時などで求める必要が出てきます。
これもフェルマーの小定理を使います。
$\ p$が素数の場合、以下の式で求める事ができます。

\begin{align*}
a^{p-1} \equiv 1 \mod{p} \\
a^{p-2} \equiv a^{-1} \mod{p} \\
\end{align*}

モジュラ逆数 - Wikipedia

平方根（square root）

$x$値から、$y$値を求める必要が出てくるので、平方根を求める方法も準備します。
素数$p$を$4$で割った余りが$3$の場合は、比較的簡単に求める事が出来ます。
secp256k1はこの条件に該当します。
まずはフェルマーの小定理から、以下を求めます。

\begin{align*}
y^{p-1} &\equiv 1 \mod{p} \\
(y^{\frac{p-1}{2}}-1)(y^{\frac{p-1}{2}}+1) &\equiv 0 \mod{p} \\
y^{\frac{p-1}{2}} &\equiv \pm 1 \mod{p} \\
\end{align*}

すこし逆説的な求め方ですが、以下のように求める事ができます。

\begin{align*}
y^2 &\equiv c \mod{p} \\
(y^2)^{\frac{p+1}{4}} &\equiv c^{\frac{p+1}{4}} \mod{p} \\
y^{\frac{p+1}{2}} &\equiv c^{\frac{p+1}{4}} \mod{p} \\
y^{\frac{p-1}{2}}\cdot y &\equiv c^{\frac{p+1}{4}} \mod{p} \\
y &\equiv \pm c^{\frac{p+1}{4}} \mod{p} \\
\end{align*}

Quadratic residue - Wikipedia

プログラミング

それでは、biginteger（任意精度整数）でelliptic curve/secp256k1（楕円曲線）の加法を実装してみます。
今回は、golangで実装しますが、どのプログラム言語でも可能だと思います。

secp256k1の素数$p$を定義します。

// p is a prime number of secp256k1.
// http://www.secg.org/sec2-v2.pdf
var p, _ = new(big.Int).SetString("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F", 16)

座標のPointを作成します。
無限遠点かどうか、複製、圧縮フォーマットを実装します。
圧縮フォーマットですが、座標$x$に対し$y$は2点あります。
剰余演算の場合、$y\mod{p}$であれば、$-y=p-y$となります。
素数は$2$以外、奇数なので$y$と$-y$は、一方が偶数で他方は奇数となります。
従って、$y$が偶数であれば0x02を、奇数であれば0x03を先頭に付与します。
33byteの固定長となります。

// Point is a coordinate of elliptic curve.
type Point struct {
	X *big.Int
	Y *big.Int
}

// Infinite returns whether it is at infinity or not.
func (point *Point) Infinite() bool {
	if point.X == nil || point.Y == nil {
		return true
	}
	return false
}

// Clone returns a copy of Point.
func (point *Point) Clone() *Point {
	clone := &Point{}
	if point.Infinite() {
		return nil
	}
	clone.X = new(big.Int).SetBytes(point.X.Bytes())
	clone.Y = new(big.Int).SetBytes(point.Y.Bytes())
	return clone
}

// Compressed returns the compressed Point.
func (point *Point) Compressed() []byte {
	if point.Infinite() {
		return nil
	}
	size := len(p.Bytes())
	bs := new(big.Int).Mod(point.X, p).Bytes()
	for len(bs) != size {
		bs = append([]byte{0x00}, bs...)
	}
	if point.Y.Bit(0) == 0 {
		bs = append([]byte{0x02}, bs...)
	} else {
		bs = append([]byte{0x03}, bs...)
	}
	return bs
}

byte配列や文字列から、Pointを生成する関数を実装します。
$x$座標から$y$座標を求める時に、前節で説明した平方根を用いる計算を使っています。
secp256k1の基点$G$を定義します。
そのオーダー$n$も定義します。

// Decode returns a Point from the bytes.
func Decode(bs []byte) (*Point, error) {
	size := len(p.Bytes())
	if len(bs) == 1+2*size {
		if bs[0] != 0x04 {
			return nil, fmt.Errorf("invalid format : %x", bs)
		}
		point := &Point{}
		point.X = new(big.Int).SetBytes(bs[1 : size+1])
		point.Y = new(big.Int).SetBytes(bs[size+1:])
		return point, nil
	}
	if len(bs) != 1+size {
		return nil, fmt.Errorf("invalid length : %x", bs)
	}
	if bs[0] != 0x02 && bs[0] != 0x03 {
		return nil, fmt.Errorf("invalid format : %x", bs)
	}
	point := &Point{}
	point.X = new(big.Int).SetBytes(bs[1:])
	// (x^3 + 7)^((p + 1) / 4)
	point.Y = new(big.Int).Exp(
		new(big.Int).Add(new(big.Int).Exp(point.X, big.NewInt(3), p), big.NewInt(7)),
		new(big.Int).Div(new(big.Int).Add(p, big.NewInt(1)), big.NewInt(4)),
		p)
	if (bs[0] != 0x02 && point.Y.Bit(0) == 0) || (bs[0] != 0x03 && point.Y.Bit(0) == 1) {
		point.Y.Sub(p, point.Y)
	}
	return point, nil
}

// DecodeString returns a Point from the hexstring.
func DecodeString(hexstring string) (*Point, error) {
	bs, err := hex.DecodeString(hexstring)
	if err != nil {
		return nil, err
	}
	return Decode(bs)
}

// G is the base point of secp256k1.
var G, _ = DecodeString("0279BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798")

// n is the order of G.
var n, _ = new(big.Int).SetString("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141", 16)

Pointの加算と倍算を実装します。
加算は、前の説で説明した計算方法で求めています。
逆数が出てくるので、前節で説明した計算で求めています。
倍算は、加算を利用しています。
※：ライブラリは射影を使っているので、計算速度はもっと速いはずです。

// Add returns the addition of Points.
func Add(P, Q *Point) *Point {
	if P.Infinite() {
		return Q.Clone()
	}
	if Q.Infinite() {
		return P.Clone()
	}
	if P.X.Cmp(Q.X) == 0 && Q.Y.Cmp(Q.Y) != 0 {
		return &Point{}
	}
	var s *big.Int
	if P.X.Cmp(Q.X) == 0 && P.Y.Cmp(Q.Y) == 0 {
		// (3xP^2) * (2 * yP)^(p - 2) mod p
		s = new(big.Int).Mod(
			new(big.Int).Mul(
				new(big.Int).Mul(new(big.Int).Mul(big.NewInt(3), P.X), P.X),
				new(big.Int).Exp(
					new(big.Int).Mul(big.NewInt(2), P.Y),
					new(big.Int).Sub(p, big.NewInt(2)),
					p)),
			p)
	} else {
		// (yP - yQ) * (xP - xQ)^(p - 2) mod p
		s = new(big.Int).Mod(
			new(big.Int).Mul(
				new(big.Int).Sub(P.Y, Q.Y),
				new(big.Int).Exp(
					new(big.Int).Sub(P.X, Q.X),
					new(big.Int).Sub(p, big.NewInt(2)),
					p)),
			p)
	}
	R := &Point{}
	// xR = s*s - (xP + xQ) mod p
	R.X = new(big.Int).Mod(new(big.Int).Sub(new(big.Int).Mul(s, s), new(big.Int).Add(P.X, Q.X)), p)
	// -yR = s*(xP - xR) - yP mod p
	R.Y = new(big.Int).Mod(new(big.Int).Sub(new(big.Int).Mul(s, new(big.Int).Sub(P.X, R.X)), P.Y), p)
	return R
}

// Mul is the multiple of Point.
func Mul(x *big.Int, P *Point) *Point {
	R := &Point{}
	for i := 0; i < x.BitLen(); i++ {
		if x.Bit(i) == 1 {
			R = Add(R, P)
		}
		P = Add(P, P)
	}
	return R
}

これで実装完了です。
正しいかどうか、Test Vectorsでチェックしてみましょう。

さいごに

Githubにサンプルコードがありますので、
詳しく見たい方はこちらをご覧ください。
また、ほとんどの言語で可能ですので一度試してみてください。