0
大学数学基礎解説
AsirでSchoofアルゴリズム実装してみた
28
0
$$$$
動機
SEAアルゴリズムの高速化に挑戦しようと思い，まずはSchoofアルゴリズムを作ってみた。
理論
省略。例えば、Lawrence C. Washington; "Elliptic Curves: Number Theory and Cryptography", CRC press, 2008 に書いてある。
実装
注意点
- コメント文の英語が大体おかしい
- Mathlogにコピペするときにインデントがうまく移っていないところがある
- 使わない関数も書いてある
- 有限体上の平方根計算を有限体上多項式の既約分解でごまかしている箇所がある
- 作動試験を十分してないのでバグはあるかも(というかたぶんある)
コード
/*** Schoof algorithm ***/
/* only if the coefficients of EC
over finite field(Fq) are over base field(Fp) */
/** sub function **/
/* improve function "var" and "vars" */
/*
Input:
F: polynomial
Output:
M: monopolynomial
*/
def varz(F) {
if (var(F) == 0) { return x; }
else { return var(F);}
}
def varsz(F) {
if (vars(F) == []) { return [x]; }
else { return vars(F);}
}
/* compute L-th "reduced" division polynomial of elliptic curve over finite field */
/*
Input:
L: integer (>=3),
degree of div poly.
EC: vector (size: 2),
[a,b] are on (Fp)^2 and coefficients of elliptic curve over Fq,
i.e., weierstrass form of elliptic curve is y^2=x^3+a*x+b.
Char: integer (prime),
characteristic of base field.
Mode: 0 or 1,
you can choose between list or polynomial as output.
Output:
if Mode = 0,
Divpoly: univariate polynomial over FF,
f_L which is L-th "reduced" division polynomial over FF,
s.t. if P=(x,y) is in E[L] and not in E[2],
then L*P=O if and only if f_L(x)=0.
if Mode = 1,
List: list (size: L+1)
list of all the division polynomials of E,
each rank of which is from 0 to L.
however, 0-th one is 0, 1-st one is 1,
2-nd one is 1.
*/
def divpoly(L, EC, Char, Mode) {
/* Caution: in the case that the caracteristic of this finite field Char
is enough small ((L^2-1)/2 >= Char), the degree of divpoly is larger
than Char, so you may have error some computation. */
if ((L^2-1)/2 >= Char) {
print("Caution: in the case that the caracteristic of
this finite field Char is enough small ((L^2-1)/2 >= Char),
the degree of divpoly is larger than Char, so you may have
error some computation.");
}
/* coefficients of EC */
A = EC[0];
B = EC[1];
ECpoly = (x^3 + A*x + B) % Char;
/* compute reduced division polynomial */
/* Case: L <= 2 */
if (L <= 2) {
error("you should not be L <= 2");
} else {}
/* Case: 2 < L <= 4 */
List_Divpoly = [
0 % Char,
1 % Char,
1 % Char,
(3*x^4+6*A*x^2+12*B*x-A^2) % Char,
(2*(x^6+5*A*x^4+20*B*x^3-5*A^2*x^2-4*A*B*x-8*B^2-A^3)) % Char
];
/*output */
if (L <= 4) {
if (Mode == 0) {
return List_Divpoly[L];
} else {
return List_Divpoly;
}
} else {}
/* Case: L >= 5 */
Vec_Divpoly = newvect(L+1,List_Divpoly);
F_sq = (16*ECpoly^2) % Char;
N = 5;
while (N <= L) {
/* Case: N = 2M+1 */
if (N % 2 == 1) {
M = (N-1)/2;
/* Case: M is odd */
if (M % 2 == 1) {
Vec_Divpoly[N] = Vec_Divpoly[M+2]*Vec_Divpoly[M]^3
- F_sq*Vec_Divpoly[M-1]*Vec_Divpoly[M+1]^3;
Vec_Divpoly[N] = Vec_Divpoly[N] % Char;
} else { /* Case: M is even */
Vec_Divpoly[N] = F_sq*Vec_Divpoly[M+2]*Vec_Divpoly[M]^3
- Vec_Divpoly[M-1]*Vec_Divpoly[M+1]^3;
Vec_Divpoly[N] = Vec_Divpoly[N] % Char;
}
} else {
/* Case: N = 2M */
M = N/2;
Vec_Divpoly[N] = (Vec_Divpoly[M+2]*Vec_Divpoly[M-1]^2
- Vec_Divpoly[M-2]*Vec_Divpoly[M+1]^2)*Vec_Divpoly[M];
Vec_Divpoly[N] = Vec_Divpoly[N] % Char;
}
N++;
}
/* Output */
if (Mode == 0) {
return Vec_Divpoly[L];
} else {
return vtol(Vec_Divpoly);
}
}
/* represent N by binary */
/*
Input:
N: integer (>=0).
Output:
BinVect: vector (length>1),
BinVect = [b0,..,bn], each bi = 0 or 1,
and N = b0*2^0 + .. + bn*2^n.
*/
def dec_to_bin (N) {
/* input need to be positive or 0 */
if (N < 0) {error("N needs to be more than 0.");} else {}
/* if input is 0, return 0 */
if (N == 0) {return vect([0]);} else {}
/* compute the digit of N and Fact_2 = max{ 2^i | N > 2^i} */
M = N;
BinDigit = 1;
Fact_2 = 1;
while (1) {
if (M == Fact_2) {
/* if N is factoring of 2 */
BinVect = newvect(BinDigit);
BinVect[BinDigit-1] = 1;
return BinVect;
} else {
/* if N is NOT factoring of 2 */
if (M < Fact_2) {
BinDigit--;
Fact_2 = Fact_2/2;
break;
} else {
/* repeat */
BinDigit++;
Fact_2 = Fact_2*2;
}
}
}
/* represent N by binary
by making M and Fact_2 smaller and smaller */
BinVect = newvect(BinDigit);
BinVect[BinDigit-1] = 1;
/* make M more smaller and find the digits which is 1 */
M = M - Fact_2;
BinDigit--;
Fact_2 = Fact_2/2;
while (M != 0) {
if (M >= Fact_2) {
/* if the digits which is 1 */
BinVect[BinDigit-1] = 1;
M = M - Fact_2;
} else {}
Fact_2 = Fact_2/2;
BinDigit--;
}
return BinVect;
}
/* compute N-th power of A by polynomial time */
/*
Input:
A: integer, polynomial
N: integer (>=0).
Output:
E: integer, polynomial
E = A^N, if both A and N are 0, error.
*/
def fast_exp (A,N) {
/* input N need to be positive or 0 */
if (N < 0) {error("input N needs to be more than 0.");} else {}
/* both input A and N are 0, error */
if (A == 0 && N == 0) {error("both input A and N are 0.");} else {}
/* if input N is 0, return 0 */
if (N == 0) {return 1;} else {}
/* if input N is 1, return A */
if (N == 1) {return A;} else {}
/* Now we represent N to binary */
BinVect = dec_to_bin(N);
/* let BinVect be [b0,b1,b2,..,bn],
compute [A,A^2,A^4,..,A^(2^n)],
and A^N = b1*A + b2*A^2 +..+ bn*A^n */
Vec_A = newvect(length(BinVect));
/* compute A^(2^0) */
Vec_A[0] = A;
Exp_A = 1;
if (BinVect[0] == 1) {Exp_A = Exp_A * Vec_A[0];}
/* compute A^(2^I) */
for (I = 1; I < length(Vec_A); I++) {
Vec_A[I] = Vec_A[I-1]^2;
if (BinVect[I] == 1) {
Exp_A = Exp_A * Vec_A[I];
}
}
return Exp_A;
}
/* compute N-th power of A mod R by polynomial time */
/*
Input:
A: integer, unipolynomial, multipolynomial
N: integer (>=0).
R: integer, unipolynomial, multipolynomial
Mode: integer (0 or 1 or 2)
if over Z, Mode = 0,
if over F[x], Mode = 1,
if over F[x1,..,xn], Mode = 2.
Output:
E: integer, unipolynomial, multipolynomial
E = (A mod R)^N,
if both A and N are 0, error.
if R = 0, error.
*/
def fast_exp_mod (A,N,R,Mode) {
/* input N need to be positive or 0 */
if (N < 0) {error("input N needs to be more than 0.");} else {}
/* both input A and N are 0, error */
if (A == 0 && N == 0) {error("both input A and N are 0.");} else {}
/* we need R is NOT 0 */
if (R == 0) {error("we need that R is NOT 0.");} else {}
/* if input N is 0, return 0 */
if (N == 0) {return 1;} else {}
/*compute B = A mod R */
if (Mode == 0) {
B = A % R;
} else {
if (Mode == 1) {
B = urem(A,R);
} else {
/* if Mode = 2 */
B = srem(A,R);
}
}
/* if input N is 1, return A */
if (N == 1) {
return B;
} else {}
/* Now we represent N to binary */
BinVect = dec_to_bin(N);
/* let BinVect be [b0,b1,b2,..,bn],
compute [A,A^2,A^4,..,A^(2^n)],
and A^N = b1*A + b2*A^2 +..+ bn*A^n */
Vec_A = newvect(length(BinVect));
/* compute A^(2^0) */
Vec_A[0] = B;
Exp_A = 1;
if (BinVect[0] == 1) {Exp_A = Exp_A * Vec_A[0];}
/* compute A^(2^I) */
for (I = 1; I < length(Vec_A); I++) {
Vec_A[I] = Vec_A[I-1]^2;
/* reduce R */
if (Mode == 0) {
Vec_A[I] = Vec_A[I] % R;
} else {
if (Mode == 1) {
Vec_A[I] = urem(Vec_A[I],R);
} else {
/* if Mode = 2 */
Vec_A[I] = srem(Vec_A[I],R);
}
}
if (BinVect[I] == 1) {
Exp_A = Exp_A * Vec_A[I];
/* reduce R */
if (Mode == 0) {
Exp_A = Exp_A % R;
} else {
if (Mode == 1) {
Exp_A = urem(Exp_A,R);
} else {
/* if Mode = 2 */
Exp_A = srem(Exp_A,R);
}
}
}
}
return Exp_A;
}
/* compute N*P,
P is a point of elliptic curve over finite field(Fq) */
/* non-used */
/*
Input:
N: integer,
P: vector (size: 3) or 0,
homogeneous coodinates of a point of EC over FF.
EC: vector (size: 2),
[a,b] are in Fp and coefficients of elliptic curve,
i.e., weierstrass form of elliptic curve is y^2=x^3+a*x+b.
Output:
Q: vector (size: 3) or 0,
homogeneous coodinates of N*P.
*/
def ecm_times_ff (N, P, EC) {
/* Case: N = 0 */
if (N == 0) {
return 0;
} else {}
/* Case: P = 0 */
if (P == 0) {
return 0;
} else {}
/* Case: N != 0 and P != 0 */
/* we can improve on this part by the represent of N */
if (N > 0) { /* Case: N > 0 */
R = P;
M = N;
} else { /* Case: N < 0 */
R = ecm_chsgn_ff(P);
M = -N;
}
/* repeat plussing R */
Q = R;
for (I = 1; I < M; I++) {
Q = ecm_add_ff(Q,R,EC);
}
return Q;
}
/* compute N*P using the binary represent of N,
P is a point of elliptic curve over finite field(Fq) */
/* non-used */
/*
Input:
N: integer,
P: vector (size: 3) or 0,
homogeneous coodinates of a point of EC over FF.
EC: vector (size: 2),
[a,b] are in Fp and coefficients of elliptic curve,
i.e., weierstrass form of elliptic curve is y^2=x^3+a*x+b.
Output:
Q: vector (size: 3) or 0,
homogeneous coodinates of N*P.
*/
def ecm_times_ff_fast (N, P, EC) {
/* Case: N = 0 */
if (N == 0) {
return 0;
} else {}
/* Case: P = 0 */
if (P == 0) {
return 0;
} else {}
/* Case: N != 0 and P != 0 */
/* we can improve on this part by the represent of N */
if (N > 0) { /* Case: N > 0 */
R = P;
M = N;
} else { /* Case: N < 0 */
R = ecm_chsgn_ff(P);
M = -N;
}
/* if input N is 1 or -1, return P or -P */
if (M == 1) {return R;} else {}
/* repeat plussing R */
/* Now we represent M to binary */
BinVect = dec_to_bin(M);
/* let BinVect be [b0,b1,b2,..,bn],
compute [R,2R,4R,..,(2*n)R],
and M*R = b1*R + b2*2R +..+ bn*(2*n)R */
Vec_R = newvect(length(BinVect));
/* compute (2^0)R */
Vec_R[0] = R;
Exp_R = 0;
if (BinVect[0] == 1) {Exp_R = ecm_add_ff(Exp_R,Vec_R[0],EC);}
/* compute (2^I)R */
for (I = 1; I < length(Vec_R); I++) {
Vec_R[I] = ecm_add_ff(Vec_R[I-1],Vec_R[I-1],EC);
if (BinVect[I] == 1) {
Exp_R = ecm_add_ff(Exp_R,Vec_R[I],EC);
}
}
return Exp_R;
}
/* extended Euclid algorithm for integer*/
/*
Input:
A: integer.
B: integer.
Output:
V: vector (size: 3),
V=(Z,X,Y) s.t. A*X+B*Y=Z=GCD(A,B).
*/
def ext_igcd (A, B) {
R1 = A;
X1 = 1;
X2 = 0;
R2 = B;
X2 = 0;
Y2 = 1;
while (R2 != 0) {
R3 = R1 % R2;
Q = idiv(R1, R2);
X3 = X1 - Q * X2;
Y3 = Y1 - Q * Y2;
R1 = R2;
R2 = R3;
X1 = X2;
X2 = X3;
Y1 = Y2;
Y2 = Y3;
}
return newvect(3,[R1, X1, Y1]);
}
/* Chinese Remainder Theory */
/*
Input:
V: vector,
V = [[a1 n1] [a2 n2] ... [am nm]],
each (ni,nj) is coprime,
we solve equations that x1 = a1 (mod n1),..,xm = am (mod nm).
Output:
S: vector (size: 2),
S = [as ns], as is the root of the equations,
ns is n1*n2*...*nm.
*/
def crt (V) {
A1 = V[0][0];
N1 = V[0][1];
for ( I = 0; I+1 < length(V); I++) {
A2 = V[I+1][0];
N2 = V[I+1][1];
EEA = ext_igcd(N1,N2);
N3 = N1*N2;
A3 = (A1*N2*EEA[2]+A2*N1*EEA[1])%N3;
A1 = A3;
N1 = N3;
}
return newvect(2,[A1%N1,N1]);
}
/* extended Euclid algorithm for unipoly over Fp */
/*
Input:
A: univariate polynomial over Fp.
B: univariate polynomial over Fp.
Char: integer (prime),
characteristic of base field.
Output:
V: vector (size: 3),
V=(Z,X,Y) s.t. A*X+B*Y=Z=GCD(A,B).
*/
def ext_ugcd (A, B, Char) {
setmod_ff(Char);
R1 = simp_ff(A);
X1 = simp_ff(1);
X2 = simp_ff(0);
R2 = simp_ff(B);
X2 = simp_ff(0);
Y2 = simp_ff(1);
while (R2 != simp_ff(0)) {
Sqr = sqr(R1,R2);
Q = Sqr[0];
R3 = Sqr[1];
X3 = X1 - Q * X2;
Y3 = Y1 - Q * Y2;
R1 = R2;
R2 = R3;
X1 = X2;
X2 = X3;
Y1 = Y2;
Y2 = Y3;
}
return newvect(3,[R1, X1, Y1]);
}
/* compute inverse of A(X) over Fp[X]/<R(X)> */
/*
Input:
Char: integer (prime),
characteristic of base field.
R: univariate polynomial over Fp.
A: univariate polynomial over Fp.
Output:
B: univariate polynomial over Fp.
if there exist the inverse of A,
then output is B such that A*B = 1.
otherwise, that is GCD(A(X),R(X)) = 1,
output is "noinverse".
*/
def inv_unipoly (A, R, Char) {
setmod_ff(Char);
/* check that there exist the inverse of A,
that is, check that GCD(A(X),R(X)) = 1 */
V = ext_ugcd(R,A,Char);
if (deg(simp_ff(V[0]),varz(simp_ff(V[0]))) > 0) {
return "noinverse";
} else {
return simp_ff(V[2])/simp_ff(V[0]);
}
}
/* compute x^N over Fp[x]/<F(X)>. */
/* non-used */
/*
Input:
Char: integer (prime),
characteristic of base field.
N: integer.
F: univariate polynomial over Fp.
Output:
R: univariate polynomial over Fp.
R is equal to x^N reduced by F(X).
*/
def red_xN_by_F (N, F, Char) {
/* set finite field */
setmod_ff(Char);
/* compute reducing */
R = urem(simp_ff(varz(F)^N), simp_ff(F));
return R;
}
/* compute Frobenius map from P = (X(x),y*Y(x)) */
/*
Input:
X: univariate polynomial,
X(x) is the polynomial which represents the x-coodinate of P.
Char: integer (prime),
characteristic of base field.
Defpoly_FF: univariate polynomial over FF,
definition polynomial of FF.
if FF is prime field, then Defpoly_FF = x.
R: univariate polynomial over Fp.
Output:
Xqn: univariate polynomial over Fp.
Frob(P) = (Xqn(x),y*Yqn(x))
*/
def frob_x(X, Char, Defpoly_FF, R) {
/* set finite field */
setmod_ff(Char);
/* compute order of field */
Order = Char^deg(Defpoly_FF,varz(Defpoly_FF));
/* reduce x^Order by R to reduce repeating */
Red_x_Order = urem(simp_ff(varz(X)^Order),simp_ff(R));
/* compute X^Order using reduced x^Order */
X_exp_Q = subst(Red_x_Order,varz(Red_x_Order),X);
X_exp_Q = urem(simp_ff(X_exp_Q),simp_ff(R));
return simp_ff(X_exp_Q);
}
/* compute Frobenius map from P = (X(x),y*Y(x)) */
/*
Input:
Y: univariate polynomial,
Y(x) is the polynomial which represents the y-coodinate of P.
EC: vector (size: 2),
[a,b] are on (Fp)^2 and coefficients of elliptic curve over Fq,
i.e., weierstrass form of elliptic curve is y^2=x^3+a*x+b.
Char: integer (prime),
characteristic of base field.
Defpoly_FF: univariate polynomial over FF,
definition polynomial of FF.
if FF is prime field, then Defpoly_FF = x.
R: univariate polynomial over Fp.
Output:
Yqn: univariate polynomial over Fp.
Frob(P) = (Xqn(x),y*Yqn(x))
*/
/* We want to rewrite this code to use reduced y^((q-1)/2)*x^q by R,
However I can't reduce y^((q-1)/2) by R */
def frob_y(Y, EC, Char, Defpoly_FF, R) {
/* set finite field */
setmod_ff(Char);
/* compute order of field */
Order = Char^deg(Defpoly_FF,varz(Defpoly_FF));
/* compute exp(y*Y,Order) = y*exp(ECpoly,(Order-1)/2)*exp(Y,Order) */
ECpoly = simp_ff(varz(Y)^3+EC[0]*varz(Y)+EC[1]);
Yp_y = fast_exp_mod(simp_ff(ECpoly),(Order-1)/2,simp_ff(R),1);
Yp_Y = fast_exp_mod(simp_ff(Y),Order,simp_ff(R),1);
Yp = urem(simp_ff(Yp_y*Yp_Y),simp_ff(R));
return simp_ff(Yp);
}
/* compute X(x), which represent (X(x),y*Y(x)) = n*(x,y) over Fp[x]/<R(x)> */
/*
Input:
N: integer (>= 1).
EC: vector (size: 2),
[a,b] are on (Fp)^2 and coefficients of elliptic curve over Fq,
i.e., weierstrass form of elliptic curve is y^2=x^3+a*x+b.
Char: integer (prime),
characteristic of base field.
R: univariate polynomial over Fp,
{the roots of R} \cap {the root of N-th divpoly} = \emptyset.
List_Divpoly: list (size: L+1)
list of all the division polynomials of E,
each rank of which is from 0 to L.
however, 0-th one is 0, 1-st one is 1,
2-nd one is 1.
Output:
X: univariate polynomial over Fp.
*/
def n_times_x (N, EC, Char, R, List_Divpoly) {
/* set finite field */
setmod_ff(Char);
/* if N = 1 */
if (N == 1) {
return simp_ff(varz(R));
} else {}
/* if N >= 2 */
ECpoly = simp_ff(varz(R)^3+EC[0]*varz(R)+EC[1]);
/* compute division polynomials */
Vec_Divpoly = ltov(List_Divpoly);
/* if N is even */
if (N % 2 == 0) {
/*
X = varz(R) - (1/4)*Vec_Divpoly[N+1]*Vec_Divpoly[N-1]
*inv_unipoly(ECpoly*(Vec_Divpoly[N]^2),R,Char);
*/
X = varz(R) - (1/4)*urem(Vec_Divpoly[N+1]*Vec_Divpoly[N-1],R)
*inv_unipoly(ECpoly*urem(Vec_Divpoly[N]^2,R),R,Char);
X = urem(X,R);
} else { /* if N is odd */
/*
X = varz(R) - 4*Vec_Divpoly[N+1]*Vec_Divpoly[N-1]
*ECpoly*inv_unipoly(Vec_Divpoly[N]^2,R,Char);
*/
X = varz(R) - 4*urem(Vec_Divpoly[N+1]*Vec_Divpoly[N-1],R)
*ECpoly*inv_unipoly(urem(Vec_Divpoly[N]^2,R),R,Char);
X = urem(X,R);
}
return simp_ff(X);
}
/* compute Y(x), which represent (X(x),y*Y(x)) = n*(x,y) over Fp[x]/<R(x)> */
/*
Input:
N: integer (>= 1).
EC: vector (size: 2),
[a,b] are on (Fp)^2 and coefficients of elliptic curve over Fq,
i.e., weierstrass form of elliptic curve is y^2=x^3+a*x+b.
Char: integer (prime),
characteristic of base field.
R: univariate polynomial over Fp.
List_Divpoly: list (size: L+1)
list of all the division polynomials of E,
each rank of which is from 0 to L.
however, 0-th one is 0, 1-st one is 1,
2-nd one is 1.
Output:
Y: univariate polynomial over Fp.
*/
def n_times_y (N, EC, Char, R, List_Divpoly) {
/* set finite field */
setmod_ff(Char);
/* if N = 1 */
if (N == 1) {
return simp_ff(1);
} else {}
/* if N >= 2 */
ECpoly = simp_ff(varz(R)^3+EC[0]*varz(R)+EC[1]);
/* compute division polynomials */
Vec_Divpoly = ltov(List_Divpoly);
/* if N is even */
if (N % 2 == 0) {
/*
Y = (1/16)*(Vec_Divpoly[N+2]*Vec_Divpoly[N-1]^2-Vec_Divpoly[N-2]*Vec_Divpoly[N+1]^2)
*inv_unipoly(ECpoly^2*Vec_Divpoly[N]^3,R,Char);
*/
Y = (1/16)*urem(Vec_Divpoly[N+2]*urem(Vec_Divpoly[N-1]^2,R)-Vec_Divpoly[N-2]*urem(Vec_Divpoly[N+1]^2,R),R)
*inv_unipoly(ECpoly^2*urem(Vec_Divpoly[N]*urem(Vec_Divpoly[N]^2,R),R),R,Char);
Y = urem(Y,R);
} else { /* if N is odd */
/*
Y = urem((urem(urem(Vec_Divpoly[N+2],R)*urem(Vec_Divpoly[N-1]^2,R),R)-urem(urem(Vec_Divpoly[N-2],R)*urem(Vec_Divpoly[N+1]^2,R),R))
*inv_unipoly(Vec_Divpoly[N]*urem(Vec_Divpoly[N]^2,R),R,Char)
,R);
*/
Y = urem(Vec_Divpoly[N+2]*urem(Vec_Divpoly[N-1]^2,R)-Vec_Divpoly[N-2]*urem(Vec_Divpoly[N+1]^2,R),R)
*inv_unipoly(Vec_Divpoly[N]*urem(Vec_Divpoly[N]^2,R),R,Char);
Y = urem(Y,R);
}
return simp_ff(Y);
}
/* compute XS(x),YS(x) which satisfy (XS(x),y*YS(x)) = (XP(x),y*YP(x)) + (XQ(x),y*YQ(x)).
in this time, P=(XP,YP) and Q=(XQ,YQ) are in E[L], L is not 2. */
/*
Input:
XP: univariate polynomial,
XP(x) is the polynomial which represents the x-coodinate of P.
YP: univariate polynomial,
YP(x) is the polynomial which represents the y-coodinate of P.
XQ: univariate polynomial,
XQ(x) is the polynomial which represents the x-coodinate of Q.
YQ: univariate polynomial,
YQ(x) is the polynomial which represents the y-coodinate of Q.
EC: vector (size: 2),
[a,b] are on (Fp)^2 and coefficients of elliptic curve over Fq,
i.e., weierstrass form of elliptic curve is y^2=x^3+a*x+b.
Char: integer (prime),
characteristic of base field.
R: univariate polynomial over Fp.
Output:
XS: univariate polynomial,
XS(x) is the polynomial which represents the x-coodinate of S = P + Q.
YS: univariate polynomial,
YS(x) is the polynomial which represents the y-coodinate of S = P + Q.
*/
def add_x_y(XP,YP,XQ,YQ,EC,Char,R) {
/* set finite field */
setmod_ff(Char);
ECpoly = simp_ff(x^3+EC[0]*x+EC[1]);
Xp = urem(simp_ff(XP),simp_ff(R));
Yp = urem(simp_ff(YP),simp_ff(R));
Xq = urem(simp_ff(XQ),simp_ff(R));
Yq = urem(simp_ff(YQ),simp_ff(R));
/* if P = Q */
if (Xp == Xq && Yp == Yq) {
Lambda = simp_ff((3*Xp^2+EC[0])*(1/2)*inv_unipoly(Yp,simp_ff(R),Char));
XS = urem(Lambda^2,simp_ff(R))*inv_unipoly(ECpoly,simp_ff(R),Char)-2*Xp;
XS = urem(simp_ff(XS),simp_ff(R));
YS = Lambda*(Xp-XS)*inv_unipoly(ECpoly,simp_ff(R),Char)-Yp;
YS = urem(simp_ff(YS),simp_ff(R));
} else {
/* if P != Q */
Lambda = simp_ff((Yq-Yp)*inv_unipoly(Xq-Xp,simp_ff(R),Char));
Lambda = urem(Lambda,simp_ff(R));
XS = urem(Lambda^2,simp_ff(R))*ECpoly-Xq-Xp;
XS = urem(simp_ff(XS),simp_ff(R));
YS = Lambda*(Xp-XS)-Yp;
YS = urem(simp_ff(YS),simp_ff(R));
}
return [simp_ff(XS),simp_ff(YS)];
}
/* compute a root of a polynomial over prime field */
/*
Input:
F: univariate polynomial over Fp.
we allow that F has no roots over Fp.
Char: integer (prime),
characteristic of base field.
Mode: integer,
0 or 1 or 2.
Output:
Root: integer or string,
if F has no roots over Fp, output is "noroot".
if F has at least one root over Fp,
Mode is 0, then return 1,
Mode is 1, then return one of the roots,
Mode is 2, then return a list which includes all roots.
*/
def uniroot_ff (F, Char, Mode) {
/* set finite field */
setmod_ff(Char);
/* check that F has at least one root over Fp */
Fpoly = simp_ff(F);
Charpoly = simp_ff(var(Fpoly)^Char-var(Fpoly));
GCD = gcd(Fpoly,Charpoly);
if (deg(GCD,var(GCD)) == 0) {
/* if F has no roots over Fp */
return "noroot";
} else {
/* if F has at least one root over Fp */
if (Mode == 0) {
return 1;
} else {
L = []; /* L includes all roots */
/* compute roots */
Fact = fctr_ff(GCD);
/* we rarely have the following cause unknown error;
internal error (SEGV)
stopped in monic_randpoly_ff at line 549 in file "/home/ayuhiss/OpenXM/lib/asir/fff"
549 S += random_ff()*V^I;
*/
for (I = 0; I < length(Fact); I++) {
if (deg(Fact[I][0],var(Fact[I][0])) == 1) {
if (Mode == 1) {
/* if Mode = 1, return a root */
return simp_ff(var(Fact[I][0])-Fact[I][0]);
} else {
/* if Mode = 2, add to the list */
L = cons(simp_ff(var(Fact[I][0])-Fact[I][0]),L);
}
} else {}
}
/* output */
if (L == []) {
return "noroot";
} else {
return L;
}
}
}
}
/* compute r,s s.t. p-1=s*2^r */
/*
Input:
P: integer (>2, prime)
Output:
L: list,
R: integer,
S: integer.
*/
def decomp_r_s (P) {
Q = P-1;
R = 0;
S = 0;
/* repeat dividing 2 */
while (Q % 2 == 0) {
Q = Q/2 ;
R++;
}
S = Q;
return [R,S];
}
/* quadratic residue test */
/*
Input:
A: integer (>=0),
P: integer (odd prime).
// Note: if the type of inputs is finitefield,
you get error at "B=B%Q" //
Output:
Test: integer (1 or -1 or 0),
output (a/p), then,
if there is a quadratic residue for a mod p, return 1,
and if there is NOT, return -1,
and if p divides a, return 0.
*/
def quad_red_test(A,P) {
/* Preparation */
/* Note: if the type of inputs is finitefield,
you get error at "B=B%Q" */
B = A;
Q = P;
if (pari(isprime,Q) == -1){error("P is not prime!");} else {}
Test = 1;
/* repeat until B is -1,0,1,2 or product of odd primes */
while (1) {
/* if B >= Q, B is replaced by B%Q */
if (B >= Q) {
B = B % Q; /* restart */
} else {
/* if B < Q */
/* if B = 0, that is B|P, return 0 */
if (B == 0) { return 0; }
else {
/* if B != 0 */
/* if B = 1,2,-1 or a odd prime */
if (B == 1 || B == 2 || B == Q-1 || (pari(isprime,B) && B > 2)) {
if (B == 1) { return Test; }
else {
if (B == 2) {
if (Q % 8 == 1 || Q % 8 == 7) { return Test; }
else { return Test*(-1); }
} else {
if (B == Q-1) { return Test*(-1); }
else {
/* if B is a odd prime, switch B and Q */
if (B % 4 == 3 && Q % 4 == 3) {
Test = (-1)*Test;
} else {}
B_dush = B;
B = Q;
Q = B_dush;
}
}
}
} else {
/* if B != 1 and 2 and -1 and a odd prime, decompose B = S*2^R */
/* and we output (B/Q) = (S/Q)*(2/Q)^R */
Lrs = decomp_r_s(B+1);
R = Lrs[0];
S = Lrs[1];
/* we use 2nd supplementary law of quadratic reciprocity*/
/* (2/Q)=(-1)^((Q^2-1)/8) */
if (Q % 8 == 1 || Q % 8 == 7) { Sup2 = 1; }
else { Sup2 = -1; }
Test = Sup2^R;
if (S == 1 || S == Q-1 || (pari(isprime,S) && S > 2)) {
/* switch S and B, and restart */
B = S;
} else {
/* this case is the most naive, compute S^((Q-1)/2) */
Test_dush = fast_exp_mod(S,(Q-1)/2,Q,0);
if (Test_dush == 1) { Test_dush = 1; }
else {
if (Test_dush+1 == Q) { Test_dush = -1; }
else{error("fail due to B^(Q-1) != 1 and -1 mod Q");}
}
return Test*Test_dush;
}
}
}
}
}
}
/* compute a square root of A in finite field
by using Tonelli-Shanks Algorithm,
computational complexity is O(log^2(p)) */
/*
Input:
A: integer.
Char: integer (prime,>=5),
characteristic of base field.
Mode: integer,
if Mode = 0, we don't quadratic test,
if Mode = 1, we do quadratic test.
Output:
X: integer or string,
if A is a quadratic residue, A = X^2 (mod Char),
which square roots is output is random,
otherwise, return "nonquad".
*/
def sqrt_ff_tonelli(A, Char, Mode) {
if (Mode == 1) {
/* check A is quadratic residue mod Char */
/* Note: if the type of inputs is finitefield,
you get error at "B=B%Q" */
Test = quad_red_test(A,Char);
if (Test == -1) {return "nonquad";}
else {
if(Test == 0) {return 0;}
else {/* if Test == 1, do continue */}
}
} else {/* we don't check A is quadratic residue mod Char */}
/* set the characteristic of finite field */
P = Char;
/* Divise into cases */
if (P % 4 == 3) {
X = fast_exp_mod(A,(P+1)/4,P,0);
return X;
} else {
if (P % 8 == 5) {
QuarRed = fast_exp_mod(A,(P-1)/4,P,0);
if (QuadRed == 1) {
X = fast_exp_mod(A,(P+3)/8,P,0);
} else {
X = fast_exp_mod(2,(P-1)/4,P,0)*fast_exp_mod(A,(P+3)/8,P,0);
}
return X;
} else {/* if P % 8 = 1, use Tonelli-Shanks Algorithm */}
}
/* Tonelli-Shanks Algorithm */
/* compute r,s s.t. p-1=s*2^r */
Lrs = decomp_r_s(Char);
R = Lrs[0];
S = Lrs[1];
/* find non-quadratic residue n mod P */
Test = 1;
N = 2;
while (Test == 1) {
Test = quad_red_test(N,P);
if (Test == -1) {
break;
} else {
if ( Test == 0) {
error("we could not find quadratic residue mod P");
} else {
/* if Test = 1, repeat */
N = pari(nextprime,N+1);
}
}
}
/* compute M = N^S */
M = fast_exp_mod(N,S,P,0);
/* compute B,U,X */
/* Step 0 */
B = fast_exp_mod(A,S,P,0);
X = fast_exp_mod(A,(S+1)/2,P,0);
U = 0;
B_pow = B;
while (B_pow != 1) {
B_pow = fast_exp_mod(B_pow,2,P,0);
U++;
}
/* Step i */
while (U != 0) {
B = (fast_exp_mod(M,2^(R-U),P,0) * B) % P;
X = (fast_exp_mod(M,2^(R-U-1),P,0) * X) % P;
U = 0;
B_pow = B;
while (B_pow != 1) {
B_pow = fast_exp_mod(B_pow,2,P,0);
U++;
}
}
return X;
}
/* compute square root of A in finite field
by finding a root of x^2-A */
/*
Input:
A: integer.
Char: integer (prime),
characteristic of base field.
Output:
Sqrt: integer or string
if A is a quadratic residue, A = Sqrt^2 (mod Char),
otherwise, return "nonquad".
*/
def sqrt_ff_quadpoly(A, Char) {
/* set finite field */
setmod_ff(Char);
/* check x^2-A has roots */
F = simp_ff(x^2 - A);
Root = uniroot_ff(F,Char,1);
if (Root == "noroot") {
return "nonquad";
} else {
return simp_ff(Root);
}
}
/** main function **/
/* count points over elliptic curve over Fq */
/*
Input:
EC: vector (size: 2),
[a,b] are on (Fp)^2 and coefficients of elliptic curve over Fq,
i.e., weierstrass form of elliptic curve is y^2=x^3+a*x+b.
Char: integer (prime),
characteristic of base field.
Defpoly_FF: univariate polynomial over FF,
definition polynomial of FF.
if FF is prime field, then Defpoly_FF = x.
Output:
N: integer ( > 0 ),
the number of points on EC over FF.
*/
def schoof(EC, Char, Defpoly_FF) {
/* set the order of finite field Fq */
Q = Char^deg(Defpoly_FF,varz(Defpoly_FF));
/* coefficients of EC */
A = EC[0]; B = EC[1];
ECpoly = (varz(Defpoly_FF)^3 + A*varz(Defpoly_FF) + B) % Char;
/* what we have to do is only finding T s.t. Q+1+T = 0 */
/* compute the range of T by Hasse's Theorem */
Range = 4*(isqrt(Q)+1);
/* compute Lm := min{prime L|\Pi_{prime \ell<=L}(\ell) > Range} */
/* make a list Prime_List = {2,3,..,Lm} */
M = 1;
L = 1;
Prime_List = [];
while (M < Range) {
L = pari(nextprime,L+1);
Prime_List = cons(L,Prime_List);
M = M*L;
}
Prime_List = reverse(Prime_List);
print("Prime_List = ",0)$print(Prime_List)$
/* compute TL = T (mod L) for each small prime L */
V = []; /* hold [TL1,L1], [TL2,L2],..,[TLm,Lm] */
/* Case: L = 2 */
L = Prime_List[0];
print("L = ",0)$print(L)$
/* check ECpoly has at least one root over Fq */
Root = uniroot_ff(ECpoly,Q,0); /* teh second input of "uniroot_ff" needs prime, but we cheat here */
if (Root == "noroot") {
/* if ECpoly has no roots in Fq */
TL = 1;
} else {
/* if ECpoly has at least one root in Fq */
TL = 0;
}
print(" TL = ",0)$print(TL)$
V = cons(newvect(2,[TL%L,L]),V);
/* compute a list of 2-th ~ (Lm+2)-th Divpoly's */
List_Divpoly = divpoly(Prime_List[length(Prime_List)-1]+2,EC,Char,1);
/* Case: L > 2 */
/* in this case, we have to find TL
by computing Frob^2(P)+Q*P+TL*Frob(P) = 0, P in E[L] */
for (J = 1; J < length(Prime_List); J++) {
/* compute a point in E[L](Fp) */
L = Prime_List[J];
print("L = ",0)$print(L)$
QL = Q % L;
Divpoly = (List_Divpoly[L]) % Char;
print("Divpoly = ",0)$print(Divpoly)$
/* check there is a point in E[L](Fp) */
Roots = uniroot_ff(Divpoly,Char,2);
print("Roots = ",0)$print(Roots)$
if (Roots == "noroot") {
Check_P = "nopoint";
} else {
Check_P = "nopoint";
for (I = 0; I < length(Roots); I++) {
Sqrt = sqrt_ff_quadpoly(subst(ECpoly,var(ECpoly),Roots[I]),Char);
//Sqrt = sqrt_ff_tonelli(subst(ECpoly,var(ECpoly),Roots[I]),Char,1);
if (Sqrt != "nonquad") {
Check_P = [Roots[I],Sqrt];
print("Check_P = ",0)$print(Check_P)$
break;
} else {}
}
}
if (Check_P != "nopoint") {
/* Case1: E[L](Fp) != \emptyset */
/* in this case, we have to find TL
s.t. (Q+1+TL)*P = 0 (P in E[L]), that is, Q+1+TL = 0 mod L. */
print("Case1: there is at least one point in E[L](Fp).")$
/* compute LHS */
LHS = (QL+1) % L;
/* find TL s.t. Q+1+TL = 0 mod L */
for (I = 0; I < L; I++) {
if ((LHS+I) % L == 0) {
TL = I; break;
} else {}
}
} else {
/* Case2: E[L](Fp) = \emptyset */
/* in this case, we have to find TL
s.t. Frob^2(P)+Q*P = -TL*Frob(P), P=(x,y) in E(Fp)\E[L]. */
print("Case2: there is no points in E[L](Fp).")$
/* compute LHS = (Frob^2)(P)+q*P */
P = [var(Divpoly),1];
Frob_P = [frob_x(P[0], Char, Defpoly_FF, Divpoly),
frob_y(P[1], EC, Char, Defpoly_FF, Divpoly)];
Frob2_P = [frob_x(Frob_P[0], Char, Defpoly_FF, Divpoly),
frob_y(Frob_P[1], EC, Char, Defpoly_FF, Divpoly)]; /* if using "subst", we cannot reduce by Divpoly */
QP = [n_times_x(QL, EC, Char, Divpoly, List_Divpoly),
n_times_y(QL, EC, Char, Divpoly, List_Divpoly)];
if (Frob2_P[0] != QP[0]) { /* over Fp */
/* Case3: (Frob^2)(P) != q*P and (Frob^2)(P) != -q*P */
print("Case3: (Frob^2)(P) != q*P and (Frob^2)(P) != -q*P.")$
LHS = add_x_y(Frob2_P[0],Frob2_P[1],QP[0],QP[1],EC,Char,Divpoly);
RHS = Frob_P;
for (I = 1; I < L; I++) {
/* if (x-cood. of LHS) = (x-cood. of RHS) */
if (LHS[0] == RHS[0]) { /* over Fp */
/* if (y-cood. of LHS) = (y-cood. of RHS) */
if (LHS[1] == RHS[1]){ /* over Fp */
TL = (-1)*I; break;
} else {
TL = I; break;
}
} else {
/* if (x-cood. of LHS) != (x-cood. of RHS) */
RHS = add_x_y(RHS[0],RHS[1],Frob_P[0],Frob_P[1],EC,Char,Divpoly);
}
if (I == L-1) {error("we could not find TL s.t. Frob^2(P)+Q*P = -TL*Frob(P).");}else{}
}
} else {
/* Case4: (Frob^2)(P) = Q*P or (Frob^2)(P) = -Q*P */
print("Case4: (Frob^2)(P) = Q*P or (Frob^2)(P) = -Q*P.")$
if (Frob2_P[1] == QP[1]) { /* over Fp */
/* Case5: (Frob^2)(P) = Q*P */
print("Case5: (Frob^2)(P) = Q*P.")$
W = sqrt_ff_tonelli(Q,L,0); /* or use "sqrt_ff_quadpoly" */
print("Now x(Frob(P)) = x(W*P). ",0)$
print("QL = ",0)$print(QL,0)$
print(", W = ",0)$print(W)$
if (W == "nonquad") {error("we could not find a square root of Char mod L.");}
/* if Frob(P) = W*P, TL = -2*W,
else if Frob(P) = -W*P, TL = 2*W. */
WPy = n_times_y(W, EC, Char, Divpoly, List_Divpoly);
if (Frob_P[1] == WPy) { /* over Fp */
TL = -2*W;
} else {
TL = 2*W;
}
} else {
/* Case6: (Frob^2)(P) = -Q*P */
print("Case6: (Frob^2)(P) = -Q*P.")$
TL = 0;
}
}
}
print(" TL = ",0)$print(TL)$
V = cons(newvect(2,[TL,L]),V);
}
/* compute TM = T mod M by CRT */
print("V = ",0)$print(V)$
V = ltov(V);
CRT = crt(V);
TM = CRT[0];
print ("CRT = ",0)$print(CRT)$
/* compute T */
if (TM <= 2*isqrt(Q)) {
T = TM;
} else {
T = TM-CRT[1];
}
print("T = ",0)$print(T,0)$print(" mod ",0)$print(M)$
/* compute #E(Fq) by Hasse's Theorem */
N = Q + 1 + T;
return N;
}
/*** executing room ***/
/* Example 1 */
/* Input */
Char = 19$ /* prime */
Defpoly_FF = x$
EC = newvect(2,[2,1])$
/* compute */
N = schoof(EC,Char,Defpoly_FF)$
/* Output */
print("\nExample 1:")$
print("Char = ",0)$print(Char)$
print("Defpoly_FF =",0)$print(Defpoly_FF)$
print("EC = ",0)$print(EC)$
print("N = ",0)$print(N)$
print("\n")$
/* Example 2 */
/* Input */
Char = 101$ /* prime */
Defpoly_FF = x$
EC = newvect(2,[3,4])$
/* compute */
N = schoof(EC,Char,Defpoly_FF)$
/* Output */
print("\nExample 2:")$
print("Char = ",0)$print(Char)$
print("Defpoly_FF =",0)$print(Defpoly_FF)$
print("EC = ",0)$print(EC)$
print("N = ",0)$print(N)$
print("\n")$
/* Example 3 */
/* Input */
Char = 1031$ /* prime */
Defpoly_FF = x$
EC = newvect(2,[0,3])$
/* compute */
N = schoof(EC,Char,Defpoly_FF)$
/* Output */
print("\nExample 3:")$
print("Char = ",0)$print(Char)$
print("Defpoly_FF =",0)$print(Defpoly_FF)$
print("EC = ",0)$print(EC)$
print("N = ",0)$print(N)$
print("\n")$
/* command in pari/GP
E = ellinit([a,b],p); N = ellcard(E);
*/
end$
実行結果
[2999] load("test.rr");
Prime_List = [2,3,5]
L = 2
TL = 1
Caution: in the case that the caracteristic of
this finite field Char is enough small ((L^2-1)/2 >= Char),
the degree of divpoly is larger than Char, so you may have
error some computation.
L = 3
Divpoly = 3*x^4+12*x^2+12*x+15
Roots = [8]
Check_P = [8,4]
Case1: there is at least one point in E[L](Fp).
TL = 1
L = 5
Divpoly = 5*x^12+10*x^10+17*x^8+5*x^7+x^6+9*x^5+12*x^4+2*x^3+5*x^2+8*x+8
Roots = noroot
Case2: there is no points in E[L](Fp).
Case3: (Frob^2)(P) != q*P and (Frob^2)(P) != -q*P.
TL = 2
V = [[ 2 5 ],[ 1 3 ],[ 1 2 ]]
CRT = [ 7 30 ]
T = 7 mod 30
Example 1:
Char = 19
Defpoly_FF =x
EC = [ 2 1 ]
N = 27
Prime_List = [2,3,5,7]
L = 2
TL = 0
L = 3
Divpoly = 3*x^4+18*x^2+48*x+92
Roots = noroot
Case2: there is no points in E[L](Fp).
Case3: (Frob^2)(P) != q*P and (Frob^2)(P) != -q*P.
TL = -1
L = 5
Divpoly = 5*x^12+85*x^10+5*x^9+65*x^8+52*x^7+79*x^6+93*x^5+28*x^4+60*x^3+53*x^2+58*x+48
Roots = noroot
Case2: there is no points in E[L](Fp).
Case4: (Frob^2)(P) = Q*P or (Frob^2)(P) = -Q*P.
Case6: (Frob^2)(P) = -Q*P.
TL = 0
L = 7
Divpoly = 7*x^24+15*x^22+20*x^21+78*x^20+70*x^19+63*x^18+47*x^17+69*x^16+53*x^15+76*x^14+74*x^13+49*x^12+49*x^11+5*x^10+76*x^9+44*x^8+22*x^7+80*x^6+45*x^5+97*x^4+65*x^3+87*x^2+77*x+94
Roots = [14,27,9]
Case2: there is no points in E[L](Fp).
Case3: (Frob^2)(P) != q*P and (Frob^2)(P) != -q*P.
TL = -3
V = [[ -3 7 ],[ 0 5 ],[ -1 3 ],[ 0 2 ]]
CRT = [ 200 210 ]
T = -10 mod 210
Example 2:
Char = 101
Defpoly_FF =x
EC = [ 3 4 ]
N = 92
Prime_List = [2,3,5,7]
L = 2
TL = 0
L = 3
Divpoly = 3*x^4+36*x
Roots = [884,0]
Check_P = [0,154]
Case1: there is at least one point in E[L](Fp).
TL = 0
L = 5
Divpoly = 5*x^12+109*x^9+933*x^6+102*x^3+915
Roots = noroot
Case2: there is no points in E[L](Fp).
Case4: (Frob^2)(P) = Q*P or (Frob^2)(P) = -Q*P.
Case6: (Frob^2)(P) = -Q*P.
TL = 0
L = 7
Divpoly = 7*x^24+491*x^21+561*x^18+807*x^15+32*x^12+133*x^9+621*x^6+385*x^3+53
Roots = noroot
Case2: there is no points in E[L](Fp).
Case4: (Frob^2)(P) = Q*P or (Frob^2)(P) = -Q*P.
Case6: (Frob^2)(P) = -Q*P.
TL = 0
V = [[ 0 7 ],[ 0 5 ],[ 0 3 ],[ 0 2 ]]
CRT = [ 0 210 ]
T = 0 mod 210
Example 3:
Char = 1031
Defpoly_FF =x
EC = [ 0 3 ]
N = 1032
今後
- 理論部分をまとめる
- SEAを作る
- SEAの改良
投稿日：2020年12月12日
この記事を高評価した人
この記事を高評価した人
高評価したユーザはいません
この記事に送られたバッジ
この記事に送られたバッジ
バッジはありません。
投稿者
投稿者
ayuhiss
0
28
コメント
コメント
他の人のコメント
コメントはありません。
読み込み中...
読み込み中