use another GLV algo.

dev
MITSUNARI Shigeo 8 years ago
parent 242e962849
commit 1d2f1160d0
  1. 136
      include/mcl/bn.hpp
  2. 125
      test/glv_test.cpp

@ -206,54 +206,46 @@ struct MapToT {
template<class Fp>
struct GLV {
typedef mcl::EcT<Fp> G1;
Fp w; // (-1 + sqrt(-3)) / 2
Fp rw; // rw = 1 / w = (-1 - sqrt(-3)) / 2
int m;
mpz_class v0, v1;
mpz_class B[2][2];
mpz_class r;
mpz_class z;
bool isNegative;
mpz_class v; // 6z^2 + 4z + 1 > 0
mpz_class c; // 2z + 1
void init(const mpz_class& r, const mpz_class& z, bool isNegative)
void init(const mpz_class& r, const mpz_class& z)
{
if (!Fp::squareRoot(w, -3)) throw cybozu::Exception("GLV:init");
w = (w - 1) / 2;
if (!Fp::squareRoot(rw, -3)) throw cybozu::Exception("GLV:init");
rw = (-1 - rw) / 2;
this->r = r;
this->z = z;
this->isNegative = isNegative;
v = 1 + z * (4 + z * 6);
c = 2 * z + 1;
m = gmp::getBitSize(r);
m = (m + fp::UnitBitSize - 1) & ~(fp::UnitBitSize - 1);// a little better size
v0 = ((6 * z * z + 4 * z + 1) << m) / r;
v1 = ((-2 * z - 1) << m) / r;
B[0][0] = 6 * z * z + 2 * z;
B[0][1] = -2 * z - 1;
B[1][0] = -2 * z - 1;
B[1][1] = -6 * z * z - 4 * z - 1;
}
/*
(p^2 mod r) (x, y) = (wx, -y)
lambda = 36z^4 - 1
lambda (x, y) = (rw x, y)
*/
void mulP2(G1& Q, const G1& P) const
void mulLambda(G1& Q, const G1& P) const
{
Fp::mul(Q.x, P.x, w);
Fp::neg(Q.y, P.y);
Fp::mul(Q.x, P.x, rw);
Q.y = P.y;
Q.z = P.z;
}
/*
s = ap^2 + b mod r
assume(s < r);
lambda = 36 z^4 - 1
x = a + b * lambda mod r
*/
void getAB(mpz_class& a, mpz_class& b, const mpz_class& s) const
void split(mpz_class& a, mpz_class& b, const mpz_class& x) const
{
assert(0 < s && s < r);
/*
s = s1 * v + s2 // s1 = s / v, s2 = s % v
= s1 * c * p^2 + s2 // vP = cp^2 P
= (s3 * v + s4) * p^2 + s2 // s3 = (s1 * c) / v, s4 = (s1 * c) % v
= (s3 * c * p^2 + s4) * p^2 + s2
= (s3 * c) * p^4 + s4 * p^2 + s2 // s5 = s3 * c, p^4 = p^2 - 1
= s5 * (p^2 - 1) + s4 * p^2 + s2
= (s4 + s5) * p^2 + (s2 - s5)
*/
mpz_class t;
mcl::gmp::divmod(a, t, s, v); // a = t / v, t = t % v
a *= c;
mcl::gmp::divmod(b, a, a, v); // b = a / v, a = a % v
b *= c;
a += b;
b = t - b;
t = (x * v0) >> m;
b = (x * v1) >> m;
a = x - (t * B[0][0] + b * B[1][0]);
b = - (t * B[0][1] + b * B[1][1]);
}
void mul(G1& Q, G1 P, mpz_class x, bool constTime = false) const
{
@ -263,72 +255,32 @@ struct GLV {
return;
}
if (x < 0) {
// G1::neg(P, P);
// x = -x;
x += r;
}
mpz_class a, b;
getAB(a, b, x);
// Q = (ap^2 + b)P
split(a, b, x);
G1 A;
mulP2(A, P);
if (a < 0) {
G1::neg(A, P);
a = -a;
} else {
A = P;
}
if (b < 0) {
G1::neg(Q, P);
b = -b;
G1::neg(P, P);
}
assert(a >= 0);
assert(b >= 0);
#if 0
G1::mul(A, A, a);
G1 B;
G1::mul(B, P, b);
G1::add(Q, A, B);
return;
#endif
#if 0 // slow
size_t nA = mcl::gmp::getBitSize(a);
size_t nB = mcl::gmp::getBitSize(b);
size_t n = std::max(nA, nB);
assert(n > 0);
G1 tbl[16];
tbl[0].clear();
tbl[1] = A;
G1::dbl(tbl[2], tbl[1]);
G1::add(tbl[3], tbl[2], tbl[1]);
for (int i = 1; i < 4; i++) {
for (int j = 0; j < 4; j++) {
G1::add(tbl[i * 4 + j], tbl[(i - 1) * 4 + j], P);
}
}
for (int i = 1; i < 16; i++) {
tbl[i].normalize();
}
if (n & 1) {
n--;
bool ai = mcl::gmp::testBit(a, n);
bool bi = mcl::gmp::testBit(b, n);
unsigned int idx = bi * 4 + ai;
Q = tbl[idx];
if (n == 0) return;
} else {
Q.clear();
}
for (int i = (int)n - 2; i >= 0; i -= 2) {
G1::dbl(Q, Q);
G1::dbl(Q, Q);
bool a0 = mcl::gmp::testBit(a, i + 0);
bool a1 = mcl::gmp::testBit(a, i + 1);
bool b0 = mcl::gmp::testBit(b, i + 0);
bool b1 = mcl::gmp::testBit(b, i + 1);
unsigned int c = b1 * 8 + b0 * 4 + a1 * 2 + a0;
if (c > 0) {
Q += tbl[c];
}
Q = P;
}
mulLambda(Q, Q);
#if 0
G1::mulBase(A, A, a);
G1::mulBase(Q, Q, b);
Q += A;
#else
A.normalize();
P.normalize();
G1 tbl[4] = { A, A, P, A + P }; // tbl[0] : dummy
Q.normalize();
G1 tbl[4] = { A, A, Q, A + Q }; // tbl[0] : dummy
tbl[3].normalize();
typedef mcl::fp::Unit Unit;
const int aN = (int)mcl::gmp::getUnitSize(a);
@ -442,7 +394,7 @@ struct ParamT {
G2::init(0, b_div_xi, mcl::ec::Proj);
G2::setOrder(r);
mapTo.init(2 * p - r);
glv.init(r, z, isNegative);
glv.init(r, z);
Fp2::pow(g[0], xi, (p - 1) / 6); // g = xi^((p-1)/6)
for (size_t i = 1; i < gN; i++) {

@ -13,6 +13,103 @@ using namespace mcl::bn256;
#define PUT(x) std::cout << #x "=" << (x) << std::endl;
/*
Skew Frobenius Map and Efficient Scalar Multiplication for Pairing-Based Cryptography
Y. Sakemi, Y. Nogami, K. Okeya, H. Kato, Y. Morikawa
*/
struct oldGLV {
Fp w; // (-1 + sqrt(-3)) / 2
mpz_class r;
mpz_class v; // 6z^2 + 4z + 1 > 0
mpz_class c; // 2z + 1
void init(const mpz_class& r, const mpz_class& z)
{
if (!Fp::squareRoot(w, -3)) throw cybozu::Exception("GLV:init");
w = (w - 1) / 2;
this->r = r;
v = 1 + z * (4 + z * 6);
c = 2 * z + 1;
}
/*
(p^2 mod r) (x, y) = (wx, -y)
*/
void mulP2(G1& Q, const G1& P) const
{
Fp::mul(Q.x, P.x, w);
Fp::neg(Q.y, P.y);
Q.z = P.z;
}
/*
x = ap^2 + b mod r
assume(x < r);
*/
void split(mpz_class& a, mpz_class& b, const mpz_class& x) const
{
assert(0 < x && x < r);
/*
x = s1 * v + s2 // s1 = x / v, s2 = x % v
= s1 * c * p^2 + s2 // vP = cp^2 P
= (s3 * v + s4) * p^2 + s2 // s3 = (s1 * c) / v, s4 = (s1 * c) % v
= (s3 * c * p^2 + s4) * p^2 + s2
= (s3 * c) * p^4 + s4 * p^2 + s2 // s5 = s3 * c, p^4 = p^2 - 1
= s5 * (p^2 - 1) + s4 * p^2 + s2
= (s4 + s5) * p^2 + (s2 - s5)
*/
mpz_class t;
mcl::gmp::divmod(a, t, x, v); // a = t / v, t = t % v
a *= c;
mcl::gmp::divmod(b, a, a, v); // b = a / v, a = a % v
b *= c;
a += b;
b = t - b;
}
template<class G1>
void mul(G1& Q, const G1& P, const mpz_class& x) const
{
G1 A, B;
mpz_class a, b;
split(a, b, x);
mulP2(A, P);
G1::mul(A, A, a);
G1::mul(B, P, b);
G1::add(Q, A, B);
}
};
template<class GLV1, class GLV2>
void compareLength(const GLV1& rhs, const GLV2& lhs)
{
cybozu::XorShift rg;
int Rc = 0;
int Lc = 0;
int eq = 0;
mpz_class R0, R1, L0, L1, x;
Fr r;
for (int i = 1; i < 1000; i++) {
r.setRand(rg);
x = r.getMpz();
rhs.split(R0, R1, x);
lhs.split(L0, L1, x);
size_t R0n = mcl::gmp::getBitSize(R0);
size_t R1n = mcl::gmp::getBitSize(R1);
size_t L0n = mcl::gmp::getBitSize(L0);
size_t L1n = mcl::gmp::getBitSize(L1);
size_t Rn = std::max(R0n, R1n);
size_t Ln = std::max(L0n, L1n);
if (Rn == Ln) {
eq++;
}
if (Rn > Ln) {
Rc++;
}
if (Rn < Ln) {
Lc++;
}
}
printf("eq=%d small is better rhs=%d, lhs=%d\n", eq, Rc, Lc);
}
void testGLV(const mcl::bn::CurveParam& cp)
{
bn384init(cp);
@ -21,8 +118,14 @@ void testGLV(const mcl::bn::CurveParam& cp)
G1 P0, P1, P2;
BN::mapToG1(P0, 1);
cybozu::XorShift rg;
oldGLV oldGlv;
oldGlv.init(BN::param.r, BN::param.z);
mcl::bn::GLV<Fp> glv;
glv.init(BN::param.r, BN::param.z, BN::param.isNegative);
glv.init(BN::param.r, BN::param.z);
compareLength(glv, oldGlv);
for (int i = 1; i < 100; i++) {
BN::mapToG1(P0, i);
Fr s;
@ -33,25 +136,11 @@ void testGLV(const mcl::bn::CurveParam& cp)
CYBOZU_TEST_EQUAL(P1, P2);
glv.mul(P2, P0, ss, true);
CYBOZU_TEST_EQUAL(P1, P2);
}
for (int i = -100; i < 100; i++) {
mpz_class ss = i;
G1::mulBase(P1, P0, ss);
glv.mul(P2, P0, ss);
CYBOZU_TEST_EQUAL(P1, P2);
glv.mul(P2, P0, ss, true);
oldGlv.mul(P2, P0, ss);
CYBOZU_TEST_EQUAL(P1, P2);
}
for (int i = -100; i < 100; i++) {
mpz_class ss = i * glv.v;
G1::mulBase(P1, P0, ss);
glv.mul(P2, P0, ss);
CYBOZU_TEST_EQUAL(P1, P2);
glv.mul(P2, P0, ss, true);
CYBOZU_TEST_EQUAL(P1, P2);
}
for (int i = -100; i < 100; i++) {
mpz_class ss = -i * glv.v + -i * glv.c;
mpz_class ss = i;
G1::mulBase(P1, P0, ss);
glv.mul(P2, P0, ss);
CYBOZU_TEST_EQUAL(P1, P2);
@ -66,7 +155,7 @@ void testGLV(const mcl::bn::CurveParam& cp)
CYBOZU_TEST_AUTO(glv)
{
testGLV(mcl::bn::CurveFp254BNb);
testGLV(mcl::bn::CurveFp382_1);
testGLV(mcl::bn::CurveFp382_2);
testGLV(mcl::bn::CurveFp254BNb);
}

Loading…
Cancel
Save