PARI GPでMontgomery curvesを用いたECDSAの実装(2)

前回作成したコードではWeierstrass標準形とMontgomery curvesを同時に実行した．今回はMontgomery curvesのみで実行してみる．

\\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]);
}
\\input and output is affine point P=(x,y)
scalar(A,B,P,k,p)={
	local(P_M,pair,kP_M,k1P_M,x_k,x_k1,y_k);
	P_M=[P[1],1];
	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],
		return([0]);
	);
	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[2],x_k,x_k1,p);
	);
	return([x_k,y_k]);
}
mdbl(A,B,P,p)={
	local(x,y);
	x=Mod(((P[1]^2-1)^2)/(4*B*(P[2]^2)),p);
	y=Mod(((2*P[1]+P[1]+A)*(3*(P[1]^2)+2*A*P[1]+1))/(2*B*P[2])-(B*(3*(((P[1]^2)+2*A*P[1]+1)^3)))/((2*B*P[2])^3)-P[2],p);
	return([x,y]);
}
madd(A,B,P1,P2,p)={
	local(x,y);
	if(P1==[0],return(P2));
	if(P2==[0],return(P1));
	if(P1[1]==P2[1],
		if(P1[2]==-P2[2],
			return([0]);,
			return(mdbl(A,B,P1,p));
		);
	);
	x=Mod((B*((P2[2]-P1[2])^2))/((P2[1]-P1[1])^2)-A-P1[1]-P2[1],p);
	y=Mod(((2*P1[1]+P2[1]+A)*(P2[2]-P1[2]))/(P2[1]-P1[1])-(B*((P2[2]-P1[2])^3))/((P2[1]-P1[1])^3)-P1[2],p);
	return([x,y]);
}
{
	\\準備
	A=84;B=1;p=30011;
	P=[Mod(7768,p),Mod(9454,p)];
	l=1249;
	print("pripare");
	print("\tbase point P=",P);
	print("\torder l=",l,"\n");
	\\鍵生成
	d_B=random(l-2)+2;
	P_B=scalar(A,B,P,d_B,p);
	print("key generation");
	print("\tsecret key d_B =",d_B);
	print("\tpublic key P_B =",P_B,"\n");
	\\署名生成
	until(u!=0&&v!=0,
		r=random(l-1)+1;
		U=scalar(A,B,P,r,p);
		m=random(p);
		u=lift(lift(U[1])%l);
		v=((m+u*d_B)/r)%l;
	);
	print("signature generation");
	print("\trandom value r=",r);
	print("\tU=",U);
	print("\tm=",m);
	print("\tsignature (u,v)=(",u,",",v,")","\n");
	\\署名検証
	d=(1/v)%l;
	V=madd(A,B,scalar(A,B,P,d*m,p),scalar(A,B,P_B,d*u,p),p);
	print("signature verification");
	print("\td=",d);
	print("\tV=",V,"\n");
	\\結果表示
	print("result");
	print("\t(u,Vx)=(",u,",",lift(V[1])%l,")");
	if(u==(lift(V[1])%l),
		print("\tOK");,
		print("\tNG");
	);
}

(18:47) gp > \r m_ecdsa.gp
pripare
        base point P=[Mod(7768, 30011), Mod(9454, 30011)]
        order l=1249
key generation
        secret key d_B =1047
        public key P_B =[Mod(24356, 30011), Mod(3467, 30011)]
signature generation
        random value r=40
        U=[Mod(28450, 30011), Mod(25680, 30011)]
        m=28429
        signature (u,v)=(972,642)
signature verification
        d=1142
        V=[Mod(28450, 30011), Mod(25680, 30011)]
result
        (u,Vx)=(972,972)
        OK