diff --git a/include/mcl/bn.h b/include/mcl/bn.h
index 4221989..42c566f 100644
--- a/include/mcl/bn.h
+++ b/include/mcl/bn.h
@@ -316,6 +316,12 @@ MCLBN_DLL_API void mclBnFp_sub(mclBnFp *z, const mclBnFp *x, const mclBnFp *y);
 MCLBN_DLL_API void mclBnFp_mul(mclBnFp *z, const mclBnFp *x, const mclBnFp *y);
 MCLBN_DLL_API void mclBnFp_div(mclBnFp *z, const mclBnFp *x, const mclBnFp *y);
 
+// y is one of square root of x
+// return 0 if success else -1
+MCLBN_DLL_API int mclBnFr_squareRoot(mclBnFr *y, const mclBnFr *x);
+MCLBN_DLL_API int mclBnFp_squareRoot(mclBnFp *y, const mclBnFp *x);
+MCLBN_DLL_API int mclBnFp2_squareRoot(mclBnFp2 *y, const mclBnFp2 *x);
+
 ////////////////////////////////////////////////
 // set zero
 MCLBN_DLL_API void mclBnG1_clear(mclBnG1 *x);
diff --git a/include/mcl/impl/bn_c_impl.hpp b/include/mcl/impl/bn_c_impl.hpp
index 291ff71..ce67455 100644
--- a/include/mcl/impl/bn_c_impl.hpp
+++ b/include/mcl/impl/bn_c_impl.hpp
@@ -269,6 +269,19 @@ void mclBnFp_div(mclBnFp *z, const mclBnFp *x, const mclBnFp *y)
 	Fp::div(*cast(z),*cast(x), *cast(y));
 }
 
+int mclBnFr_squareRoot(mclBnFr *y, const mclBnFr *x)
+{
+	return Fr::squareRoot(*cast(y), *cast(x)) ? 0 : -1;
+}
+int mclBnFp_squareRoot(mclBnFp *y, const mclBnFp *x)
+{
+	return Fp::squareRoot(*cast(y), *cast(x)) ? 0 : -1;
+}
+int mclBnFp2_squareRoot(mclBnFp2 *y, const mclBnFp2 *x)
+{
+	return Fp2::squareRoot(*cast(y), *cast(x)) ? 0 : -1;
+}
+
 ////////////////////////////////////////////////
 // set zero
 void mclBnG1_clear(mclBnG1 *x)
diff --git a/test/bn_c_test.hpp b/test/bn_c_test.hpp
index 654ab90..62d7871 100644
--- a/test/bn_c_test.hpp
+++ b/test/bn_c_test.hpp
@@ -760,6 +760,75 @@ CYBOZU_TEST_AUTO(Fp2)
 	CYBOZU_TEST_ASSERT(mclBnFp2_isEqual(&x1, &x2));
 }
 
+CYBOZU_TEST_AUTO(squareRootFr)
+{
+	mclBnFr x, y, y2;
+	for (int i = 0; i < 10; i++) {
+		mclBnFr_setInt(&x, i * i);
+		CYBOZU_TEST_EQUAL(mclBnFr_squareRoot(&y, &x), 0);
+		mclBnFr_sqr(&y2, &y);
+		CYBOZU_TEST_EQUAL(mclBnFr_isEqual(&x, &y2), 1);
+	}
+	char buf[128];
+	mclBnFr_setInt(&x, -1);
+	CYBOZU_TEST_ASSERT(mclBnFr_serialize(buf, sizeof(buf), &x) > 0);
+	int mod8 = (buf[0] + 1) & 7;
+	/*
+		(2)
+		(p) = (-1)^((p^2-1)/8) = 1 if and only if there is x s.t. x^2 = 2 mod p
+	*/
+	bool hasSquareRoot = (((mod8 * mod8 - 1) / 8) & 1) == 0;
+	printf("Fr:hasSquareRoot=%d\n", hasSquareRoot);
+	mclBnFr_setInt(&x, 2);
+	CYBOZU_TEST_EQUAL(mclBnFr_squareRoot(&y, &x), hasSquareRoot ? 0 : -1);
+	if (hasSquareRoot) {
+		mclBnFr_sqr(&y2, &y);
+		CYBOZU_TEST_EQUAL(mclBnFr_isEqual(&x, &y2), 1);
+	}
+}
+
+CYBOZU_TEST_AUTO(squareRootFp)
+{
+	mclBnFp x, y, y2;
+	for (int i = 0; i < 10; i++) {
+		mclBnFp_setInt(&x, i * i);
+		CYBOZU_TEST_EQUAL(mclBnFp_squareRoot(&y, &x), 0);
+		mclBnFp_sqr(&y2, &y);
+		CYBOZU_TEST_EQUAL(mclBnFp_isEqual(&x, &y2), 1);
+	}
+	char buf[128];
+	mclBnFp_setInt(&x, -1);
+	CYBOZU_TEST_ASSERT(mclBnFp_serialize(buf, sizeof(buf), &x) > 0);
+	int mod8 = (buf[0] + 1) & 7;
+	/*
+		(2)
+		(p) = (-1)^((p^2-1)/8) = 1 if and only if there is x s.t. x^2 = 2 mod p
+	*/
+	bool hasSquareRoot = (((mod8 * mod8 - 1) / 8) & 1) == 0;
+	printf("Fp:hasSquareRoot=%d\n", hasSquareRoot);
+	mclBnFp_setInt(&x, 2);
+	CYBOZU_TEST_EQUAL(mclBnFp_squareRoot(&y, &x), hasSquareRoot ? 0 : -1);
+	if (hasSquareRoot) {
+		mclBnFp_sqr(&y2, &y);
+		CYBOZU_TEST_EQUAL(mclBnFp_isEqual(&x, &y2), 1);
+	}
+}
+
+#if 0
+CYBOZU_TEST_AUTO(squareRootFp2)
+{
+	mclBnFp2 x, y, y2;
+	for (int i = 0; i < 10; i++) {
+		mclBnFp_setByCSPRNG(&x.d[0]);
+		mclBnFp_setByCSPRNG(&x.d[1]);
+		mclBnFp2_sqr(&x, &x);
+		CYBOZU_TEST_EQUAL(mclBnFp2_squareRoot(&y, &x), 0);
+		mclBnFp2_sqr(&y2, &y);
+		CYBOZU_TEST_EQUAL(mclBnFp2_isEqual(&x, &y2), 1);
+	}
+}
+#endif
+
 CYBOZU_TEST_AUTO(mapToG1)
 {
 	mclBnFp x;