[mathjax]今回はPARI GPでMontgomery curvesの加算，2倍算の実装を行う．
もともとPARI GPにはWeierstrass標準形の関数しか用意されていないため，曲線や点を相互に変換することにより動作確認を行いつつ実装を行う．Weierstass，Montgomery間の曲線と点の変換はWikipediaに記載されている手順により行うことができる．Montgomery curvesの加算，2倍算についてもWikipediaに記載されているものを実装する．
ただし，加算については\(mP+nP\)のとき\(m>n\)かつ\(m-n=1\)として簡略化して記載する．なぜなら，Montgomery ladderを使用する際にはm-n=1が常に成り立つため，実用上は問題ないためである．
1.加算の実装
以下の式に従いコードを書く．
$$ X_{m+n}=Z_{m-n}((X_m-Z_m)(X_n+Z_n)+(X_m+Z_m)(X_n-Z_n))^2 $$
$$ Z_{m+n}=X_{m-n}((X_m-Z_m)(X_n+Z_n)-(X_m+Z_m)(X_n-Z_n))^2 $$
add(X1,Z1,Xm,Zm,Xn,Zn,p)={
local(a,b,c,d,Xmn,Zmn);
a=Xm+Zm;
b=Xm-Zm;
c=Xn+Zn;
d=Xn-Zn;
Xmn=Mod(Z1*((b*c+a*d)^2),p);
Zmn=Mod(X1*((b*c-a*d)^2),p);
return([Xmn,Zmn]);
}
2.2倍算の実装
以下の式に従いコードを書く．
$$4X_nZ_n=(X_n+Z_n)^2-(X_n-Z_n)^2$$
$$X_{2n}=(X_n+Z_n)^2(X_n-Z_n)^2$$
$$Z_{2n}=4X_nZ_n((X_n-Z_n)^2+((A+2)/4)(4X_nZ_n))$$
dbl(A,X,Z,p)={
local(tmp1,tmp2,tmp3,X2,Z2);
tmp1=(X+Z)^2;
tmp2=(X-Z)^2;
tmp3=tmp1-tmp2;
X2=Mod(tmp1*tmp2,p);
Z2=Mod(tmp3*(tmp2+((tmp3*(A+2))/4)),p);
return([X2,Z2]);
}
3.y座標の復元
動作確認のために必要なので，y座標の復元を行う．式はHandbook of Elliptic and Hyperelliptic Curve Cryptographyに記載してある以下の式を使う．
$$y_n = \frac{(x_1 x_n +1)(x_1 + x_n +2A)-2A-(x_1-x_n)^2x_{n+1}}{2By_1}$$
コードとしては以下のようになる．
recover(A,B,x1,y1,xn,xn1,p)={
local(up,down);
up=(x1*xn+1)*(x1+xn+2*A)-2*A-((x1-xn)^2)*xn1;
down=2*B*y1;
return(Mod(up/down,p));
}
4.動作確認
4.1.曲線の決定と変換
今回は例としてMontgomery curvesの1つである\(y^2=x^3+84x^2+x\)を使って動作確認を行う．最初に\(y^2=x^3+84x^2+x\)をWeierstrass標準形にする．
アルゴリズムより\(A=84,B=1\)となるため，
$$\begin{array}{rcl}a&=&\frac{3-A^2}{3B^2}\\&=&-2351\end{array}$$
$$\begin{array}{rcl}b&=&\frac{2A^3-9A}{27B^3}\\&=&43876\end{array}$$
となり，Montgomery curvesの\(y^2=x^3+84x^2+x\)はWeierstrass標準形では\(y^2=x^3-2351x+43876\)となる．
4.2.ベースポイント\(P_W\)の生成と\(2P_W\)の確認
次に，PARI GPの機能を用いてWeierstrass標準形でランダムな点\(P_W\)を生成し，\(P_W,2P_W\)を確認する．
なお，\(p=251\)とし，有限体\(F_{251}\)上の楕円曲線とする．
(17:41) gp > e=ellinit([0,0,0,-2351,43876],251) %88 = [Mod(0, 251), Mod(0, 251), Mod(0, 251), Mod(159, 251), Mod(202, 251), Mod(0, 251), Mod(67, 251), Mod(55, 251), Mod(70, 251), Mod(149, 251), Mod(168, 251), Mod(133, 251), Mod(170, 251), Vecsmall([3]), [251, [244, 215, [6, 0, 0, 0]]], [0, 0, 0, 0]] (17:46) gp > Pw=random(e) %89 = [Mod(201, 251), Mod(28, 251)] (17:46) gp > ellmul(e,Pw,2) %90 = [Mod(50, 251), Mod(154, 251)]
\(P_W=(201,28),2P_W=(50,154)\)となった．
4.3.\(P_W\)をMontgomery curvesの点\(P_M\)に変換
Weierstrass標準形の点\(P_W=(t,v)\)をMontgomery curveの点\(P_M=(x,y)\)に変換する．
\(s=1,\alpha=28\)となるため，
$$\begin{array}{rcl}(t,v)&=>&(x,y)\\&=&(s(t-\alpha),sv)\\&=&(t-28,v)\end{array}$$
$$(201,28)=>(173,28)$$
よって，\(P_W=(201,28)\)は\(P_M=(173,28)\)と変換できる．
4.4.\(P_M\)を射影座標に変換
\(P_M=(x,y)=(173,28)\)を射影座標に変換すると，\(P_M=(X_1,Y_1,Z_1)=(173,28,1)\)となる．さらにy座標は演算に必要ないためこれを簡略化して，以降は\(P_M=(X_1,Z_1)=(173,1)\)とする．
4.5.\(2P_M,3P_M\)の確認
先ほど作成したコードをつかって，\(2P_M=dbl(P_M),3P_M=add(2P_M,P_M)\)を求める．なお，\(3P_M\)を求めるのは\(2P_M\)のy座標を復元するためである．
(14:51) gp > dbl(84,173,1,251) %4 = [Mod(218, 251), Mod(124, 251)] (14:52) gp > add(173,1,218,124,173,1,251) %5 = [Mod(93, 251), Mod(219, 251)]
よって，
$$P_M=(173,1)\\2P_M=(218,124)\\3P_M=(93,219)$$
となることがわかる．
4.6.\(2P_M=(x_2,y_2)\)の\(y_2\)の復元
\(2P_M\)のy座標を復元する．その準備として，各点を射影座標からアフィン座標へと変換する．変換は\((x,y)=(X/Z,Y/Z)\)として行う．なお，\(2P_M,3P_M\)のy座標はYがわからないため計算不能で，そこは?としている．
$$\begin{array}{rcl}P_M&=&(X_1,Z_1)&=&(173,1)&,&(x_1,y_1)&=&(173,28)\\2P_M&=&(X_2,Z_2)&=&(218,124)&,&(x_2,y_2)&=&(22,?)\\3P_M&=&(X_3,Z_3)&=&(93,219)&,&(x_3,y_3)&=&(52,?)\end{array}$$
\(x_1,y_1,x_2,x_3\)が揃ったので\(y_2\)の復元を行う．
(15:22) gp > recover(84,1,173,28,22,52,251) %12 = Mod(154, 251)
その結果，\(2P_M=(x_2,y_2)=(22,154)\)となることがわかった．
4.7.\(2P_M\)が\(2P_W\)と等しくなることの確認
\(2P_M=(x_2,y_2)=(22,154)=>(22+28,154)=(50,154)=2P_W\)となるため，点の2倍算は正しく実行できたことが分かる．また，y座標を復元するために点の加算を行ったため，点の加算についても正しく実装できたと考える．
5.今回のソースコード
add(X1,Z1,Xm,Zm,Xn,Zn,p)={
local(a,b,c,d,Xmn,Zmn);
a=Xm+Zm;
b=Xm-Zm;
c=Xn+Zn;
d=Xn-Zn;
Xmn=Mod(Z1*((b*c+a*d)^2),p);
Zmn=Mod(X1*((b*c-a*d)^2),p);
return([Xmn,Zmn]);
}
dbl(A,X,Z,p)={
local(tmp1,tmp2,tmp3,X2,Z2);
tmp1=(X+Z)^2;
tmp2=(X-Z)^2;
tmp3=tmp1-tmp2;
X2=Mod(tmp1*tmp2,p);
Z2=Mod(tmp3*(tmp2+((tmp3*(A+2))/4)),p);
return([X2,Z2]);
}
recover(A,B,x1,y1,xn,xn1,p)={
local(up,down);
up=(x1*xn+1)*(x1+xn+2*A)-2*A-((x1-xn)^2)*xn1;
down=2*B*y1;
return(Mod(up/down,p));
}