improve squareRoot by removing legendre

Makefile

@@ -28,7 +28,7 @@ ifeq ($(MCL_STATIC_CODE),1)
 endif
 ifeq ($(CPU),x86-64)
   MCL_USE_XBYAK?=1
-  TEST_SRC+=mont_fp_test.cpp sq_test.cpp
+  TEST_SRC+=mont_fp_test.cpp #sq_test.cpp
   ifeq ($(MCL_USE_XBYAK),1)
     TEST_SRC+=fp_generator_test.cpp
   endif

include/mcl/fp.hpp

@@ -207,15 +207,7 @@ class FpT : public fp::Serializable<FpT<tag, maxBitSize>,
 	}
 	static inline bool squareRoot(FpT& y, const FpT& x)
 	{
-		if (isMont()) return op_.sq.get(y, x);
-		mpz_class mx, my;
-		bool b = false;
-		x.getMpz(&b, mx);
-		if (!b) return false;
-		b = op_.sq.get(my, mx);
-		if (!b) return false;
-		y.setMpz(&b, my);
-		return b;
+		return op_.sq.get(y, x);
 	}
 	FpT() {}
 	FpT(const FpT& x)

include/mcl/gmp_util.hpp

@@ -784,6 +784,61 @@ class SquareRoot {
 		}
 		return false;
 	}
+	/*
+		solve x^2 = a in Fp
+	*/
+	template<class Fp>
+	bool getCandidate(Fp& x, const Fp& a) const
+	{
+		assert(Fp::getOp().mp == p);
+		if (a.isZero() || a.isOne()) {
+			x = a;
+			return true;
+		}
+		if (r == 1) {
+			// (p + 1) / 4 = (q + 1) / 2
+			Fp::pow(x, a, q_add_1_div_2);
+			return true;
+		}
+		Fp c, d;
+		{
+			bool b;
+			c.setMpz(&b, s);
+			assert(b);
+		}
+		int e = r;
+		Fp::pow(d, a, q);
+		Fp::pow(x, a, q_add_1_div_2); // destroy a if &x == &a
+		Fp dd;
+		Fp b;
+		while (!d.isOne()) {
+			int i = 1;
+			Fp::sqr(dd, d);
+			while (!dd.isOne()) {
+				Fp::sqr(dd, dd);
+				i++;
+				if (i >= e) return false;
+			}
+			assert(e > i);
+			int t = e - i - 1;
+			const int tMax = 30; // int32_t max
+			if (t < tMax) {
+				b = 1 << t;
+			} else {
+				b = 1 << tMax;
+				t -= tMax;
+				for (int j = 0; j < t; j++) {
+					b += b;
+				}
+			}
+			Fp::pow(b, c, b);
+			x *= b;
+			Fp::sqr(c, b);
+			d *= c;
+			e = i;
+		}
+		return true;
+	}
 public:
 	SquareRoot() { clear(); }
 	bool isPrecomputed() const { return isPrecomputed_; }
@@ -838,100 +893,18 @@ class SquareRoot {
 		q_add_1_div_2 = (q + 1) / 2;
 		*pb = true;
 	}
-	/*
-		solve x^2 = a mod p
-	*/
-	bool get(mpz_class& x, const mpz_class& a) const
-	{
-		if (!isPrime) {
-			return false;
-		}
-		if (a == 0) {
-			x = 0;
-			return true;
-		}
-		if (gmp::legendre(a, p) < 0) return false;
-		if (r == 1) {
-			// (p + 1) / 4 = (q + 1) / 2
-			gmp::powMod(x, a, q_add_1_div_2, p);
-			return true;
-		}
-		mpz_class c = s, d;
-		int e = r;
-		gmp::powMod(d, a, q, p);
-		gmp::powMod(x, a, q_add_1_div_2, p); // destroy a if &x == &a
-		mpz_class dd;
-		mpz_class b;
-		while (d != 1) {
-			int i = 1;
-			dd = d * d; dd %= p;
-			while (dd != 1) {
-				dd *= dd; dd %= p;
-				i++;
-			}
-			b = 1;
-			b <<= e - i - 1;
-			gmp::powMod(b, c, b, p);
-			x *= b; x %= p;
-			c = b * b; c %= p;
-			d *= c; d %= p;
-			e = i;
-		}
-		return true;
-	}
-	/*
-		solve x^2 = a in Fp
-	*/
-	template<class Fp>
-	bool get(Fp& x, const Fp& a) const
+	template<class T>
+	bool get(T& x, const T& a) const
 	{
-		assert(Fp::getOp().mp == p);
-		if (a == 0) {
-			x = 0;
-			return true;
-		}
-		{
-			bool b;
-			mpz_class aa;
-			a.getMpz(&b, aa);
-			assert(b);
-			if (gmp::legendre(aa, p) < 0) return false;
-		}
-		if (r == 1) {
-			// (p + 1) / 4 = (q + 1) / 2
-			Fp::pow(x, a, q_add_1_div_2);
-			return true;
-		}
-		Fp c, d;
-		{
-			bool b;
-			c.setMpz(&b, s);
-			assert(b);
-		}
-		int e = r;
-		Fp::pow(d, a, q);
-		Fp::pow(x, a, q_add_1_div_2); // destroy a if &x == &a
-		Fp dd;
-		Fp b;
-		while (!d.isOne()) {
-			int i = 1;
-			Fp::sqr(dd, d);
-			while (!dd.isOne()) {
-				dd *= dd;
-				i++;
-			}
-			b = 1;
-//			b <<= e - i - 1;
-			for (int j = 0; j < e - i - 1; j++) {
-				b += b;
+		T t, t2;
+		if (getCandidate(t, a)) {
+			T::sqr(t2, t);
+			if (t2 == a) {
+				x = t;
+				return true;
 			}
-			Fp::pow(b, c, b);
-			x *= b;
-			Fp::sqr(c, b);
-			d *= c;
-			e = i;
 		}
-		return true;
+		return false;
 	}
 	bool operator==(const SquareRoot& rhs) const
 	{

include/mcl/vint.hpp

@@ -696,6 +696,9 @@ class Vint {
 	// logical left shift (copy sign)
 	static void shl(Vint& y, const Vint& x, size_t shiftBit)
 	{
+if (shiftBit > MCL_MAX_BIT_SIZE*2) {
+	printf("shiftBit=%zd\n", shiftBit);
+}
 		assert(shiftBit <= MCL_MAX_BIT_SIZE * 2); // many be too big
 		size_t xn = x.size();
 		size_t yn = xn + (shiftBit + UnitBitSize - 1) / UnitBitSize;

test/bench.hpp

@@ -88,6 +88,20 @@ void invVecBench(const char *msg)
 	CYBOZU_BENCH_C(msg, C, mcl::invVec, x, x, n);
 }
 
+template<class T>
+void sqrBench(const T& x, const char *msg)
+{
+	T y = x * x, z;
+	CYBOZU_BENCH_C((std::string(msg) + "::squareRoot T").c_str(), 10000, T::squareRoot, z, y);
+	for (int i = 0; i < 100; i++) {
+		if (!T::squareRoot(z, y)) {
+			break;
+		}
+		y += T::one();
+	}
+	CYBOZU_BENCH_C((std::string(msg) + "::squareRoot F").c_str(), 10000, T::squareRoot, z, y);
+}
+
 void testBench(const G1& P, const G2& Q)
 {
 #ifndef NDEBUG
@@ -151,6 +165,7 @@ void testBench(const G1& P, const G2& Q)
 	CYBOZU_BENCH_C("Fp::sqr       ", C3, Fp::sqr, x, x);
 	CYBOZU_BENCH_C("Fp::inv       ", C3, invAdd, x, x, y);
 	CYBOZU_BENCH_C("Fp::pow       ", C3, Fp::pow, x, x, y);
+	sqrBench(x, "Fp");
 	invVecBench<Fp>("Fp:invVec");
 	invVecBench<Fr>("Fr:invVec");
 	{
@@ -168,6 +183,7 @@ void testBench(const G1& P, const G2& Q)
 		CYBOZU_BENCH_C("Fr::sqr       ", C3, Fr::sqr, a, a);
 		CYBOZU_BENCH_C("Fr::inv       ", C3, invAdd, a, a, b);
 		CYBOZU_BENCH_C("Fr::pow       ", C3, Fr::pow, a, a, b);
+		sqrBench(a, "Fr");
 	}
 	Fp2 xx, yy;
 	xx.a = x;

test/ec_test.cpp

@@ -385,19 +385,37 @@ struct Test {
 	}
 	void squareRoot() const
 	{
+		Ec P(Fp(para.gx), Fp(para.gy));
+
+		for (int i = 0; i < 100; i++) {
+			Ec::dbl(P, P);
+			P.normalize();
+			Fp x = P.x;
+			Fp y = P.y;
+			Fp yy;
+			CYBOZU_TEST_ASSERT(Ec::getYfromX(yy, x, y.isOdd()));
+			CYBOZU_TEST_EQUAL(yy, y);
+			Fp::neg(y, y);
+			yy.clear();
+			CYBOZU_TEST_ASSERT(Ec::getYfromX(yy, x, y.isOdd()));
+			CYBOZU_TEST_EQUAL(yy, y);
+			yy += P.y;
+			CYBOZU_TEST_ASSERT(yy.isZero());
+		}
+
 		Fp x(para.gx);
-		Fp y(para.gy);
-		bool odd = y.isOdd();
-		Fp yy;
-		bool b = Ec::getYfromX(yy, x, odd);
-		CYBOZU_TEST_ASSERT(b);
-		CYBOZU_TEST_EQUAL(yy, y);
-		Fp::neg(y, y);
-		odd = y.isOdd();
-		yy.clear();
-		b = Ec::getYfromX(yy, x, odd);
-		CYBOZU_TEST_ASSERT(b);
-		CYBOZU_TEST_EQUAL(yy, y);
+		for (int i = 0; i < 100; i++) {
+			mpz_class mx = x.getMpz();
+			int ret = mcl::gmp::legendre(mx, Fp::getOp().mp);
+			Fp y;
+			if (Fp::squareRoot(y, x)) {
+				CYBOZU_TEST_EQUAL(y*y, x);
+				CYBOZU_TEST_EQUAL(ret, 1);
+			} else {
+				CYBOZU_TEST_EQUAL(ret, -1);
+			}
+			x += 1;
+		}
 	}
 	void mul_fp() const
 	{

test/sq_test.cpp

@@ -11,10 +11,15 @@ CYBOZU_TEST_AUTO(sqrt)
 		sq.set(p);
 		for (mpz_class a = 0; a < p; a++) {
 			mpz_class x;
-			if (sq.get(x, a)) {
+			bool b1 = sq.get(x, a);
+			mpz_class x2;
+			bool b2 = sq.get2(x2, a);
+			CYBOZU_TEST_EQUAL(b1, b2);
+			if (b1) {
 				mpz_class y;
 				y = (x * x) % p;
 				CYBOZU_TEST_EQUAL(a, y);
+				CYBOZU_TEST_EQUAL(x, x2);
 			}
 		}
 	}