PARI GPにMontgomery ladderを実装する

[mathjax]前回の続き．今回はMontgomery ladder．

今回はMontgomery ladderを実装する．アルゴリズムはWikipediaにも載っている．見れば分かるように，仕組みとしてはバイナリ法と同じ程度に簡素なものとなっている．しかし，バイナリ法と比較してサイドチャネル攻撃に強いという特性を持つ．バイナリ法ではk[i]の値が1か0かで加算の有無が決まり，これを外部から観察するとkが復元されるという弱点が会った．一方でMontgomery ladderではk[i]の値に関わらず加算と2倍算を毎回行うので，サイドチャネル攻撃に耐性を持つ．
Montgomery ladderのコードは以下のようになる．

ladder(A,P,k,p)={
	local(R0,R1,i,bin);
	bin=binary(k);
	R0=[0];
	R1=P;
	for(i=1,length(bin),
		if(bin[i]==0,
			R1=add(P,R1,R0,p);
			R0=dbl(A,R0,p);,
			R0=add(P,R1,R0,p);
			R1=dbl(A,R1,p);
		);
	);
	return([R1,R0]);
}

 
前回のコードを修正しつつMontgomery ladderのコードを足し，テスト用のコードも付けた全体のコードを以下に示す．
修正点としては，addとdblの引数をX,Zという形式から点Pの形式に変更した．また，無限遠点に関わる分岐を追加した．recoverについても\(P_{n+1}\)が無限遠点の場合の処理についてのコードを追加した．
 

\\m-n==1
add(P1,Pm,Pn,p)={
	local(a,b,c,d,Xmn,Zmn);
	if(Pm==0||Pm[2]==0,return(Pn));
	if(Pn==0||Pn[2]==0,return(Pm));
	a=Pm[1]+Pm[2];
	b=Pm[1]-Pm[2];
	c=Pn[1]+Pn[2];
	d=Pn[1]-Pn[2];
	Xmn=Mod(P1[2]*((b*c+a*d)^2),p);
	Zmn=Mod(P1[1]*((b*c-a*d)^2),p);
	return([Xmn,Zmn]);
}
dbl(A,P,p)={
	local(tmp1,tmp2,tmp3,X2,Z2);
	if(P==0||P[2]==0,return([0]));
	tmp1=(P[1]+P[2])^2;
	tmp2=(P[1]-P[2])^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,res);
	if(xn1==-1,
		return(-y1);
	);
	up=(x1*xn+1)*(x1+xn+2*A)-2*A-((x1-xn)^2)*xn1;
	down=2*B*y1;
	if(down==0,
		res=0;,
		res=Mod(up/down,p);
	);
	return(res);
}
ladder(A,P,k,p)={
	local(R0,R1,i,bin);
	bin=binary(k);
	R0=[0];
	R1=P;
	for(i=1,length(bin),
		if(bin[i]==0,
			R1=add(P,R1,R0,p);
			R0=dbl(A,R0,p);,
			R0=add(P,R1,R0,p);
			R1=dbl(A,R1,p);
		);
	);
	return([R1,R0]);
}
{
while(1,
	A=84;B=1;
	p=251;
	e=ellinit([0,0,0,-2351,43876],p);
	P_W=random(e);
	P_M=[P_W[1]-28,1];
	k=random(p);
	kP_W=ellmul(e,P_W,k);
	pair=ladder(A,P_M,k,p);
	kP_M=pair[2];k1P_M=pair[1];
	if(kP_M!=[0]&&kP_M[2]==0,kP_M=[0]);
	if(k1P_M!=[0]&&k1P_M[2]==0,k1P_M=[0]);
	if(kP_M==[0]&&kP_M==kP_W,
		/*print("GOOD!!");*/next();
	);
	if(kP_M!=[0],
		x_k=Mod(kP_M[1]/kP_M[2],p);
		if(k1P_M==[0]||k1P_M[2]==0,
			x_k1=-1;,
			x_k1=Mod(k1P_M[1]/k1P_M[2],p);
		);
		y_k=recover(A,B,P_M[1],P_W[2],x_k,x_k1,p);
		if(kP_W==[x_k+28,y_k],
			/*print("GOOD!!");*/,
			print("bad...");
			print("recover(",A,"  ",B,"  ",P_M[1],"  ",P_W[2],"  ",x_k,"  ",x_k1,"  ",p,")");
			print("P_W ",P_W);
			print("P_M ",P_M);
			print("k ",k);
			print("kP_W ",kP_W);
			print("pair ",pair	);
			print("kP_M ",kP_M);
			print("k1P_M ",k1P_M);
			print("x_k ",x_k);
			print("y_k ",y_k);
			print(kP_W,"==",[x_k+28,y_k]);
		);
	);
);
}