refactor mapto

update-fork
MITSUNARI Shigeo 5 years ago
parent 7e6ab212fb
commit a0dcc763f8
  1. 224
      include/mcl/ec.hpp
  2. 108
      include/mcl/mapto_wb19.hpp

@ -51,6 +51,11 @@ enum ModeCoeffA {
namespace local {
/*
elliptic class E must have
member variables of type Fp x, y, z
static member a_, b_, specialA_
*/
// x is negative <=> x < half(:=(p+1)/2) <=> a = 1
template<class F>
bool get_a_flag(const F& x)
@ -68,61 +73,69 @@ bool get_a_flag(const mcl::Fp2T<F>& x)
} // mcl::ec::local
template<class F>
void normalizeJacobi(F& x, F& y, F& z)
template<class E>
void normalizeJacobi(E& P)
{
if (z.isZero()) return;
F::inv(z, z);
typedef typename E::Fp F;
if (P.z.isZero()) return;
F::inv(P.z, P.z);
F rz2;
F::sqr(rz2, z);
x *= rz2;
y *= rz2;
y *= z;
z = 1;
F::sqr(rz2, P.z);
P.x *= rz2;
P.y *= rz2;
P.y *= P.z;
P.z = 1;
}
// (x/z^2, y/z^3)
template<class F>
bool isEqualJacobi(const F& x1, const F& y1, const F& z1, const F& x2, const F& y2, const F& z2)
template<class E>
bool isEqualJacobi(const E& P1, const E& P2)
{
typedef typename E::Fp F;
F s1, s2, t1, t2;
F::sqr(s1, z1);
F::sqr(s2, z2);
F::mul(t1, x1, s2);
F::mul(t2, x2, s1);
F::sqr(s1, P1.z);
F::sqr(s2, P2.z);
F::mul(t1, P1.x, s2);
F::mul(t2, P2.x, s1);
if (t1 != t2) return false;
F::mul(t1, y1, s2);
F::mul(t2, y2, s1);
t1 *= z2;
t2 *= z1;
F::mul(t1, P1.y, s2);
F::mul(t2, P2.y, s1);
t1 *= P2.z;
t2 *= P1.z;
return t1 == t2;
}
// Y^2 == X(X^2 + aZ^4) + bZ^6
template<class F>
bool isValidJacobi(const F& x, const F& y, const F& z, const F& a, const F& b)
template<class E>
bool isValidJacobi(const E& P)
{
typedef typename E::Fp F;
F y2, x2, z2, z4, t;
F::sqr(x2, x);
F::sqr(y2, y);
F::sqr(z2, z);
F::sqr(x2, P.x);
F::sqr(y2, P.y);
F::sqr(z2, P.z);
F::sqr(z4, z2);
F::mul(t, z4, a);
F::mul(t, z4, E::a_);
t += x2;
t *= x;
t *= P.x;
z4 *= z2;
z4 *= b;
z4 *= E::b_;
t += z4;
return y2 == t;
}
/*
M > S + A
a = 0 2M + 5S + 14A
a = -3 2M + 7S + 15A
generic 3M + 7S + 15A
M == S
a = 0 3M + 4S + 13A
a = -3 3M + 6S + 14A
generic 4M + 6S + 14A
*/
template<class E>
void dblJacobi(E& R, const E& P, int specialA, const typename E::Fp& a)
void dblJacobi(E& R, const E& P)
{
typedef typename E::Fp F;
if (P.isZero()) {
@ -133,13 +146,19 @@ void dblJacobi(E& R, const E& P, int specialA, const typename E::Fp& a)
F x2, y2, xy, t;
F::sqr(x2, P.x);
F::sqr(y2, P.y);
F::add(xy, P.x, y2);
F::sqr(y2, y2);
F::sqr(xy, xy);
xy -= x2;
xy -= y2;
xy += xy;
switch (specialA) {
if (sizeof(F) <= 32) { // M == S
F::mul(xy, P.x, y2);
xy += xy;
F::sqr(y2, y2);
} else { // M > S + A
F::add(xy, P.x, y2);
F::sqr(y2, y2);
F::sqr(xy, xy);
xy -= x2;
xy -= y2;
}
xy += xy; // 4xy^2
switch (E::specialA_) {
case Zero:
F::add(t, x2, x2);
x2 += t;
@ -158,11 +177,11 @@ void dblJacobi(E& R, const E& P, int specialA, const typename E::Fp& a)
case GenericA:
default:
if (isPzOne) {
t = a;
t = E::a_;
} else {
F::sqr(t, P.z);
F::sqr(t, t);
t *= a;
t *= E::a_;
}
t += x2;
x2 += x2;
@ -187,9 +206,10 @@ void dblJacobi(E& R, const E& P, int specialA, const typename E::Fp& a)
}
// 7M + 4S + 7A
template<class E, class F>
void addJacobi(E& R, const E& P, const E& Q, int specialA, const F& a)
template<class E>
void addJacobi(E& R, const E& P, const E& Q)
{
typedef typename E::Fp F;
if (P.isZero()) { R = Q; return; }
if (Q.isZero()) { R = P; return; }
bool isPzOne = P.z.isOne();
@ -230,7 +250,7 @@ void addJacobi(E& R, const E& P, const E& Q, int specialA, const F& a)
r -= S1;
if (H.isZero()) {
if (r.isZero()) {
ec::dblJacobi(R, P, specialA, a);
ec::dblJacobi(R, P);
} else {
R.clear();
}
@ -263,43 +283,46 @@ void addJacobi(E& R, const E& P, const E& Q, int specialA, const F& a)
F::sub(R.y, U1, H3);
}
template<class F>
void normalizeProj(F& x, F& y, F& z)
template<class E>
void normalizeProj(E& P)
{
if (z.isZero()) return;
F::inv(z, z);
x *= z;
y *= z;
z = 1;
typedef typename E::Fp F;
if (P.z.isZero()) return;
F::inv(P.z, P.z);
P.x *= P.z;
P.y *= P.z;
P.z = 1;
}
// (Y^2 - bZ^2)Z = X(X^2 + aZ^2)
template<class F>
bool isValidProj(const F& x, const F& y, const F& z, const F& a, const F& b)
template<class E>
bool isValidProj(const E& P)
{
typedef typename E::Fp F;
F y2, x2, z2, t;
F::sqr(x2, x);
F::sqr(y2, y);
F::sqr(z2, z);
F::mul(t, a, z2);
F::sqr(x2, P.x);
F::sqr(y2, P.y);
F::sqr(z2, P.z);
F::mul(t, E::a_, z2);
t += x2;
t *= x;
z2 *= b;
t *= P.x;
z2 *= E::b_;
y2 -= z2;
y2 *= z;
y2 *= P.z;
return y2 == t;
}
// (x/z, y/z)
template<class F>
bool isEqualProj(const F& x1, const F& y1, const F& z1, const F& x2, const F& y2, const F& z2)
template<class E>
bool isEqualProj(const E& P1, const E& P2)
{
typedef typename E::Fp F;
F t1, t2;
F::mul(t1, x1, z2);
F::mul(t2, x2, z1);
F::mul(t1, P1.x, P2.z);
F::mul(t2, P2.x, P1.z);
if (t1 != t2) return false;
F::mul(t1, y1, z2);
F::mul(t2, y2, z1);
F::mul(t1, P1.y, P2.z);
F::mul(t2, P2.y, P1.z);
return t1 == t2;
}
@ -309,16 +332,17 @@ bool isEqualProj(const F& x1, const F& y1, const F& z1, const F& x2, const F& y2
sqr| 4| 5| 5
add| 11|12|12
*/
template<class E, class F>
void dblProj(E& R, const E& P, int specialA, const F& a)
template<class E>
void dblProj(E& R, const E& P)
{
typedef typename E::Fp F;
if (P.isZero()) {
R.clear();
return;
}
const bool isPzOne = P.z.isOne();
F w, t, h;
switch (specialA) {
switch (E::specialA_) {
case Zero:
F::sqr(w, P.x);
F::add(t, w, w);
@ -338,10 +362,10 @@ void dblProj(E& R, const E& P, int specialA, const F& a)
case GenericA:
default:
if (isPzOne) {
w = a;
w = E::a_;
} else {
F::sqr(w, P.z);
w *= a;
w *= E::a_;
}
F::sqr(t, P.x);
w += t;
@ -379,9 +403,10 @@ void dblProj(E& R, const E& P, int specialA, const F& a)
sqr| 2
add| 7
*/
template<class E, class F>
void addProj(E& R, const E& P, const E& Q, int specialA, const F& a)
template<class E>
void addProj(E& R, const E& P, const E& Q)
{
typedef typename E::Fp F;
if (P.isZero()) { R = Q; return; }
if (Q.isZero()) { R = P; return; }
bool isPzOne = P.z.isOne();
@ -404,7 +429,7 @@ void addProj(E& R, const E& P, const E& Q, int specialA, const F& a)
v -= r;
if (v.isZero()) {
if (A == PyQz) {
dblProj(R, P, specialA, a);
dblProj(R, P);
} else {
R.clear();
}
@ -442,22 +467,25 @@ void addProj(E& R, const E& P, const E& Q, int specialA, const F& a)
}
// y^2 == (x^2 + a)x + b
template<class F>
bool isValidAffine(const F& x, const F& y, const F& a, const F& b)
template<class E>
bool isValidAffine(const E& P)
{
typedef typename E::Fp F;
assert(!P.z.isZero());
F y2, t;
F::sqr(y2, y);
F::sqr(t, x);
t += a;
t *= x;
t += b;
F::sqr(y2, P.y);
F::sqr(t, P.x);
t += E::a_;
t *= P.x;
t += E::b_;
return y2 == t;
}
// y^2 = x^3 + ax + b
template<class E, class F>
static inline void dblAffine(E& R, const E& P, const F& a)
template<class E>
static inline void dblAffine(E& R, const E& P)
{
typedef typename E::Fp F;
if (P.isZero()) {
R.clear();
return;
@ -470,7 +498,7 @@ static inline void dblAffine(E& R, const E& P, const F& a)
F::sqr(t, P.x);
F::add(s, t, t);
t += s;
t += a;
t += E::a_;
F::add(s, P.y, P.y);
t /= s;
F::sqr(s, t);
@ -484,16 +512,17 @@ static inline void dblAffine(E& R, const E& P, const F& a)
R.z = 1;
}
template<class E, class F>
void addAffine(E& R, const E& P, const E& Q, const F& a)
template<class E>
void addAffine(E& R, const E& P, const E& Q)
{
typedef typename E::Fp F;
if (P.isZero()) { R = Q; return; }
if (Q.isZero()) { R = P; return; }
F t;
F::sub(t, Q.x, P.x);
if (t.isZero()) {
if (P.y == Q.y) {
dblAffine(R, P, a);
dblAffine(R, P);
} else {
R.clear();
}
@ -552,17 +581,17 @@ public:
private:
bool isValidAffine() const
{
return ec::isValidAffine(x, y, a_, b_);
return ec::isValidAffine(*this);
}
public:
void normalize()
{
switch (mode_) {
case ec::Jacobi:
ec::normalizeJacobi(x, y, z);
ec::normalizeJacobi(*this);
break;
case ec::Proj:
ec::normalizeProj(x, y, z);
ec::normalizeProj(*this);
break;
}
}
@ -629,10 +658,10 @@ public:
{
switch (mode_) {
case ec::Jacobi:
if (!ec::isValidJacobi(x, y, z, a_, b_)) return false;
if (!ec::isValidJacobi(*this)) return false;
break;
case ec::Proj:
if (!ec::isValidProj(x, y, z, a_, b_)) return false;
if (!ec::isValidProj(*this)) return false;
break;
case ec::Affine:
if (z.isZero()) return true;
@ -664,26 +693,27 @@ public:
{
switch (mode_) {
case ec::Jacobi:
ec::dblJacobi(R, P, specialA_, a_);
ec::dblJacobi(R, P);
break;
case ec::Proj:
ec::dblProj(R, P, specialA_, a_);
ec::dblProj(R, P);
break;
case ec::Affine:
ec::dblAffine(R, P, a_);
ec::dblAffine(R, P);
break;
}
}
static inline void add(EcT& R, const EcT& P, const EcT& Q) {
static inline void add(EcT& R, const EcT& P, const EcT& Q)
{
switch (mode_) {
case ec::Jacobi:
ec::addJacobi(R, P, Q, specialA_, a_);
ec::addJacobi(R, P, Q);
break;
case ec::Proj:
ec::addProj(R, P, Q, specialA_, a_);
ec::addProj(R, P, Q);
break;
case ec::Affine:
ec::addAffine(R, P, Q, a_);
ec::addAffine(R, P, Q);
break;
}
}
@ -1075,9 +1105,9 @@ public:
{
switch (mode_) {
case ec::Jacobi:
return ec::isEqualJacobi(x, y, z, rhs.x, rhs.y, rhs.z);
return ec::isEqualJacobi(*this, rhs);
case ec::Proj:
return ec::isEqualProj(x, y, z, rhs.x, rhs.y, rhs.z);
return ec::isEqualProj(*this, rhs);
case ec::Affine:
default:
return x == rhs.x && y == rhs.y && z == rhs.z;

@ -43,6 +43,9 @@ template<class F>
struct PointT {
typedef F Fp;
F x, y, z;
static F a_;
static F b_;
static int specialA_;
bool isZero() const
{
return z.isZero();
@ -53,15 +56,24 @@ struct PointT {
y.clear();
z.clear();
}
#if 0
bool isEqual(const PointT<F>& rhs) const
{
return ec::isEqualJacobi(*this, rhs);
}
#endif
};
template<class F> F PointT<F>::a_;
template<class F> F PointT<F>::b_;
template<class F> int PointT<F>::specialA_;
} // mcl::local
template<class Fp, class Fp2, class G2>
struct MapToG2_WB19 {
typedef local::PointT<Fp2> Point;
Fp2 xi;
Fp2 Ell2p_a;
Fp2 Ell2p_b;
Fp half;
mpz_class sqrtConst; // (p^2 - 9) / 16
Fp2 root4[4];
@ -71,7 +83,6 @@ struct MapToG2_WB19 {
Fp2 ynum[4];
Fp2 yden[4];
int draftVersion_;
typedef local::PointT<Fp2> Point;
void setDraftVersion(int version)
{
draftVersion_ = version;
@ -174,16 +185,12 @@ struct MapToG2_WB19 {
void dbl(Point& Q, const Point& P) const
{
dblT(Q, P);
}
void dbl(G2& Q, const G2& P) const
{
dblT(Q, P);
// G2::dbl(Q, P);
// ec::dblJacobi(Q, P);
}
// P is on y^2 = x^3 + Ell2p_a x + Ell2p_b
bool isValidPoint(const Point& P) const
{
return ec::isValidJacobi(P.x, P.y, P.z, Ell2p_a, Ell2p_b);
return ec::isValidJacobi(P);
}
bool isValidPoint(const G2& P) const
{
@ -194,10 +201,11 @@ struct MapToG2_WB19 {
bool b;
xi.a = -2;
xi.b = -1;
Ell2p_a.a = 0;
Ell2p_a.b = 240;
Ell2p_b.a = 1012;
Ell2p_b.b = 1012;
Point::a_.a = 0;
Point::a_.b = 240;
Point::b_.a = 1012;
Point::b_.b = 1012;
Point::specialA_ = ec::GenericA;
half = -1;
half /= 2;
sqrtConst = Fp::getOp().mp;
@ -386,11 +394,11 @@ struct MapToG2_WB19 {
// (t^2 * xi)^2 + (t^2 * xi)
Fp2::add(deno, t2xi2, t2xi);
Fp2::add(nume, deno, 1);
nume *= Ell2p_b;
nume *= Point::b_;
if (deno.isZero()) {
Fp2::mul(deno, Ell2p_a, xi);
Fp2::mul(deno, Point::a_, xi);
} else {
deno *= -Ell2p_a;
deno *= -Point::a_;
}
Fp2 u, v;
{
@ -398,8 +406,8 @@ struct MapToG2_WB19 {
Fp2::sqr(deno2, deno);
Fp2::mul(v, deno2, deno);
Fp2::mul(u, Ell2p_b, v);
Fp2::mul(tmp, Ell2p_a, nume);
Fp2::mul(u, Point::b_, v);
Fp2::mul(tmp, Point::a_, nume);
tmp *= deno2;
u += tmp;
Fp2::sqr(tmp, nume);
@ -506,20 +514,20 @@ struct MapToG2_WB19 {
den += den2;
Fp2 x0_num, x0_den;
Fp2::add(x0_num, den, 1);
x0_num *= Ell2p_b;
x0_num *= Point::b_;
if (den.isZero()) {
Fp2::mul(x0_den, Ell2p_a, xi);
Fp2::mul(x0_den, Point::a_, xi);
} else {
Fp2::mul(x0_den, -Ell2p_a, den);
Fp2::mul(x0_den, -Point::a_, den);
}
Fp2 x0_den2, x0_den3, gx0_den, gx0_num;
Fp2::sqr(x0_den2, x0_den);
Fp2::mul(x0_den3, x0_den2, x0_den);
gx0_den = x0_den3;
Fp2::mul(gx0_num, Ell2p_b, gx0_den);
Fp2::mul(gx0_num, Point::b_, gx0_den);
Fp2 tmp, tmp1, tmp2;
Fp2::mul(tmp, Ell2p_a, x0_num);
Fp2::mul(tmp, Point::a_, x0_num);
tmp *= x0_den2;
gx0_num += tmp;
Fp2::sqr(tmp, x0_num);
@ -581,22 +589,22 @@ struct MapToG2_WB19 {
{
G2 t[16];
t[0] = P;
dbl(t[1], t[0]);
add(t[4], t[1], t[0]);
add(t[2], t[4], t[1]);
add(t[3], t[2], t[1]);
add(t[11], t[3], t[1]);
add(t[9], t[11], t[1]);
add(t[10], t[9], t[1]);
add(t[5], t[10], t[1]);
add(t[7], t[5], t[1]);
add(t[15], t[7], t[1]);
add(t[13], t[15], t[1]);
add(t[6], t[13], t[1]);
add(t[14], t[6], t[1]);
add(t[12], t[14], t[1]);
add(t[8], t[12], t[1]);
dbl(t[1], t[6]);
G2::dbl(t[1], t[0]);
G2::add(t[4], t[1], t[0]);
G2::add(t[2], t[4], t[1]);
G2::add(t[3], t[2], t[1]);
G2::add(t[11], t[3], t[1]);
G2::add(t[9], t[11], t[1]);
G2::add(t[10], t[9], t[1]);
G2::add(t[5], t[10], t[1]);
G2::add(t[7], t[5], t[1]);
G2::add(t[15], t[7], t[1]);
G2::add(t[13], t[15], t[1]);
G2::add(t[6], t[13], t[1]);
G2::add(t[14], t[6], t[1]);
G2::add(t[12], t[14], t[1]);
G2::add(t[8], t[12], t[1]);
G2::dbl(t[1], t[6]);
const struct {
uint32_t n;
@ -618,21 +626,21 @@ struct MapToG2_WB19 {
};
for (size_t j = 0; j < CYBOZU_NUM_OF_ARRAY(tbl); j++) {
const uint32_t n = tbl[j].n;
for (size_t i = 0; i < n; i++) dbl(t[1], t[1]);
add(t[1], t[1], t[tbl[j].idx]);
for (size_t i = 0; i < n; i++) G2::dbl(t[1], t[1]);
G2::add(t[1], t[1], t[tbl[j].idx]);
}
for (size_t i = 0; i < 5; i++) dbl(t[1], t[1]);
add(out, t[1], t[2]);
for (size_t i = 0; i < 5; i++) G2::dbl(t[1], t[1]);
G2::add(out, t[1], t[2]);
}
void mx_chain(G2& Q, const G2& P) const
{
G2 T;
dbl(T, P);
G2::dbl(T, P);
const size_t tbl[] = { 2, 3, 9, 32, 16 };
for (size_t i = 0; i < CYBOZU_NUM_OF_ARRAY(tbl); i++) {
add(T, T, P);
G2::add(T, T, P);
for (size_t j = 0; j < tbl[i]; j++) {
dbl(T, T);
G2::dbl(T, T);
}
}
Q = T;
@ -644,12 +652,12 @@ struct MapToG2_WB19 {
#else
G2 work, work2;
h2_chain(work, P);
dbl(work2, work);
add(work2, work, work2);
G2::dbl(work2, work);
G2::add(work2, work, work2);
mx_chain(work, work2);
mx_chain(work, work);
neg(work2, work2);
add(Q, work, work2);
G2::neg(work2, work2);
G2::add(Q, work, work2);
#endif
}
template<class T>

Loading…
Cancel
Save