add fp2

include/mcl/fp.hpp

@@ -79,6 +79,12 @@ public:
 		op_.fp_add = fp_addW;
 		op_.fp_sub = fp_subW;
 		op_.fp_mul = fp_mulW;
+		// Fp2
+		op_.fp2_add = fp2_addW;
+		op_.fp2_sub = fp2_subW;
+		op_.fp2_mul = fp2_mulW;
+		op_.fp2_inv = fp2_invW;
+		op_.fp2_sqr = fp2_sqrW;
 /*
 	priority : MCL_USE_XBYAK > MCL_USE_LLVM > none
 	Xbyak > llvm_opt > llvm > gmp
@@ -453,6 +459,106 @@ private:
 	{
 		op_.mont(y, x, x, op_.p, op_.rp);
 	}
+	/*
+		default Fp2 operator
+		Fp2 = Fp[u]/(u^2 + 1)
+	*/
+	static inline void fp2_addW(Unit *z, const Unit *x, const Unit *y)
+	{
+		const fp::void3u fp_add = op_.fp_add;
+		const size_t N = op_.N;
+		fp_add(z, x, y);
+		fp_add(z + N, x + N, y + N);
+	}
+	static inline void fp2_subW(Unit *z, const Unit *x, const Unit *y)
+	{
+		const fp::void3u fp_sub = op_.fp_sub;
+		const size_t N = op_.N;
+		fp_sub(z, x, y);
+		fp_sub(z + N, x + N, y + N);
+	}
+	static inline void fp2_negW(Unit *y, const Unit *x)
+	{
+		const fp::void2u fp_neg = op_.fp_neg;
+		const size_t N = op_.N;
+		fp_neg(y, x);
+		fp_neg(y + N, x + N);
+	}
+	/*
+		x = a + bu, y = c + du, u^2 = -1
+		z = xy = (a + bu)(c + du) = (ac - bd) + (ad + bc)u
+		ad+bc = (a + b)(c + d) - ac - bd
+	*/
+	static inline void fp2_mulW(Unit *z, const Unit *x, const Unit *y)
+	{
+		const fp::void3u fp_add = op_.fp_add;
+		const fp::void3u fp_sub = op_.fp_sub;
+		const fp::void3u fp_mul = op_.fp_mul;
+		const size_t N = op_.N;
+		const Unit *a = x;
+		const Unit *b = x + N;
+		const Unit *c = y;
+		const Unit *d = y + N;
+		Unit t1[maxSize];
+		Unit t2[maxSize];
+		Unit ac[maxSize];
+		Unit bd[maxSize];
+		fp_add(t1, a, b);
+		fp_add(t2, c, d);
+		fp_mul(t1, t1, t2); // (a + b)(c + d)
+		fp_mul(ac, a, c);
+		fp_mul(bd, b, d);
+		fp_sub(z, ac, bd); // ac - bd
+		fp_sub(z, z, bd);
+		fp_sub(z + N, t1, ac);
+		fp_sub(z + N, z + N, bd);
+	}
+	/*
+		x = a + bu, u^2 = -1
+		y = x^2 = (a + bu)^2 = (a^2 - b^2) + 2abu
+	*/
+	static inline void fp2_sqrW(Unit *y, const Unit *x)
+	{
+		const fp::void3u fp_add = op_.fp_add;
+		const fp::void3u fp_sub = op_.fp_sub;
+		const fp::void2u fp_sqr = op_.fp_sqr;
+		const fp::void3u fp_mul = op_.fp_mul;
+		const size_t N = op_.N;
+		const Unit *a = x;
+		const Unit *b = x + N;
+		Unit aa[maxSize];
+		Unit bb[maxSize];
+		Unit t[maxSize];
+		fp_sqr(aa, a);
+		fp_sqr(bb, b);
+		fp_mul(t, a, b);
+		fp_sub(y, aa, bb); // a^2 - b^2
+		fp_add(y + N, t, t); // 2ab
+	}
+	/*
+		x = a + bu
+		1 / x = (a - bu) / (a^2 + b^2)
+	*/
+	static inline void fp2_invW(Unit *y, const Unit *x)
+	{
+		const fp::void3u fp_add = op_.fp_add;
+		const fp::void2u fp_sqr = op_.fp_sqr;
+		const fp::void3u fp_mul = op_.fp_mul;
+		const fp::void2uOp fp_invOp = op_.fp_invOp;
+		const fp::void2u fp_neg = op_.fp_neg;
+		const size_t N = op_.N;
+		const Unit *a = x;
+		const Unit *b = x + N;
+		Unit aa[maxSize];
+		Unit bb[maxSize];
+		fp_sqr(aa, a);
+		fp_sqr(bb, b);
+		fp_add(aa, aa, bb);
+		fp_invOp(aa, aa, op_); // aa = 1 / (a^2 + b^2)
+		fp_mul(y, y, aa);
+		fp_mul(y + N, y + N, aa);
+		fp_neg(y + N, y + N);
+	}
 };
 
 template<class tag, size_t maxBitSize> fp::Op FpT<tag, maxBitSize>::op_;

include/mcl/op.hpp

@@ -94,6 +94,20 @@ struct Op {
 	void3u fp_mulPreP;
 	void3u fp_modP;
 	FpGenerator *fg;
+
+	/*
+		for Fp2 = F[u] / (u^2 + 1)
+		x = a + bu
+	*/
+	int xi_c; // xi = u + xi_c
+	void3u fp2_add;
+	void3u fp2_sub;
+	void3u fp2_mul;
+	void2u fp2_neg;
+	void2u fp2_inv;
+	void2u fp2_sqr;
+	void2u fp2_mul_xi;
+
 	Op()
 		: N(0), bitSize(0)
 		, fp_isZero(0), fp_clear(0), fp_copy(0)
@@ -102,6 +116,9 @@ struct Op {
 		, rp(0), mont(0)
 		, fp_negP(0), fp_sqrPreP(0), fp_invOp(0), fp_addP(0), fp_subP(0), fp_mulPreP(0), fp_modP(0)
 		, fg(createFpGenerator())
+		, xi_c(0)
+		, fp2_add(0), fp2_sub(0), fp2_mul(0), fp2_neg(0)
+		, fp2_sqr(0), fp2_mul_xi(0)
 	{
 	}
 	~Op()
@@ -124,6 +141,7 @@ struct Op {
 		fp_mul(y, x, R2);
 	}
 	void init(const std::string& mstr, int base, size_t maxBitSize, Mode mode);
+	void initFp2(int xi_c);
 	static FpGenerator* createFpGenerator();
 	static void destroyFpGenerator(FpGenerator *fg);
 private:
@@ -131,24 +149,4 @@ private:
 	void operator=(const Op&);
 };
 
-/*
-	for Fp2 = F[u] / (u^2 + 1)
-	x = a + bu
-*/
-struct Op2 {
-	Op *op;
-	int xi_c; // xi = u + xi_c
-	void3u add;
-	void3u sub;
-	void3u mul;
-	void2u neg;
-	void2u sqr;
-	void2u mul_xi;
-	Op2()
-		: op(0), xi_c(0), add(0), sub(0), mul(0), neg(0), sqr(0), mul_xi(0)
-	{
-	}
-	void init(Op *op, int xi_c);
-};
-
 } } // mcl::fp

src/fp.cpp

@@ -323,12 +323,6 @@ void Op::init(const std::string& mstr, int base, size_t maxBitSize, Mode mode)
 	sq.set(mp);
 }
 
-void Op2::init(Op *op, int xi_c)
-{
-	if (op->N * UnitBitSize != 256) throw cybozu::Exception("Op2:init:not support size") << op->N;
-	this->op = op;
-}
-
 void arrayToStr(std::string& str, const Unit *x, size_t n, int base, bool withPrefix)
 {
 	switch (base) {
@@ -431,5 +425,11 @@ int64_t getInt64(bool *pb, fp::Block& b, const fp::Op& op)
 	return 0;
 }
 
+void Op::initFp2(int xi_c)
+{
+	if (N * UnitBitSize != 256) throw cybozu::Exception("Op2:init:not support size") << N;
+	this->xi_c = xi_c;
+}
+
 } } // mcl::fp