有限域上多项式乘法的快速实现

最新推荐文章于 2023-10-22 21:04:32 发布
最新推荐文章于 2023-10-22 21:04:32 发布
阅读量1.5k 收藏 6
点赞数 2
分类专栏： 代码 文章标签： 密码学 数学 抽象代数 c++ 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_44885334/article/details/123804999
版权

问题

多项式环 R = G F [ p d ] [ x ] / ( x n − 1 ) R=GF[p^d][x]/(x^n-1) R=GF[pd][x]/(xn−1)，其中有限域 G F ( p d ) ≅ Z p [ x ] / ( m ( x ) ) GF(p^d) \cong Z_p[x]/(m(x)) GF(pd)≅Zp​[x]/(m(x))， m ( x ) m(x) m(x)是 Z p Z_p Zp​上 d d d次不可约多项式。对于任意 f , g ∈ R f,g \in R f,g∈R，如何快速计算 f ⋅ g f\cdot g f⋅g？

NTL

代码
#include <iostream>
#include "tools.h"

#include <NTL/ZZ_p.h> // integers mod p
#include <NTL/ZZ_pX.h> // polynomials over ZZ_p
#include <NTL/ZZ_pE.h> // ring/field extension of ZZ_p
#include <NTL/ZZ_pEX.h> // polynomials over ZZ_pE
#include <NTL/ZZ_pXFactoring.h>
#include <NTL/ZZ_pEXFactoring.h>

using namespace std;
using namespace NTL;

#pragma comment(lib, "NTL")

int main()
{
	ZZ p(2); //初始化为2

	//群Z_p
	ZZ_p::init(p);

	//随机生成Z_p[x]中的d次不可约多项式
	int d = 149;
	ZZ_pX m;
	//BuildIrred(m, d);
	SetCoeff(m, d);
	SetCoeff(m, 10);
	SetCoeff(m, 9);
	SetCoeff(m, 7);
	SetCoeff(m, 0);

	//域GF(p^d) = Z_p[x]/m(x)
	ZZ_pE::init(m);

	//GF(p^d)[x]中的多项式
	ZZ_pEX f, g, h, P;
	int n = 53;

	//P(x)=x^n-1
	SetCoeff(P, 0, -1);
	SetCoeff(P, n);

	// f(x) = x^8 - 1
	SetCoeff(f, 8); //将 x^8 系数置为 1
	SetCoeff(f, 0, -1); //将 x^0 系数置为 -1

	//随机生成n长多项式g(x)
	random(g, n-1);

	// 环上多项式的运算
	h = g * g;

	cout << "p = " << p << endl;
	cout << "d = " << d << endl;
	cout << "n = " << n << endl;
	cout << "m(x) = " << m << endl;
	cout << "P(x) = " << P << endl;

	//cout << "f = " << f << endl;
	//cout << "g = " << g << endl;
	//cout << "h = " << h << endl;

	int loop = 10;

	printf("%d rounds, NTL mul: ", loop);
	Loop(10, f*g);

	printf("%d rounds, NTL mod: ", loop);
	Loop(10, h%P);


	return 0;
}
结果
p = 2
d = 149
n = 53
m(x) = [1 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
P(x) = [[1] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [] [1]]
10 rounds, NTL mul: 5.20162 s
10 rounds, NTL mod: 17.3201 s

NTT+Barrett

代码
IncompleteNTT
#include "IncompleteNTT.h"


int32 brv(int32 b, int32 l)
{
	int32 bb = 0;

	for (int32 i = 0; i < l; i++)
	{
		bb <<= 1;
		bb |= (b & 1);
		b >>= 1;
	}

	return bb;
}

//int64 exgcd(int64* x, int64* y, int64 a, int64 b)
//{
//	if (b == 0)
//	{
//		*x = 1;
//		*y = 0;
//		return a;
//	}
//	int64 ret = exgcd(x, y, b, a%b);
//	int64 tmp = *x;
//	*x = *y;
//	*y = tmp - (a / b) * (*y);
//	return ret;
//}

//int32 deg(Poly& f)
//{
//	int32 d = f.size() - 1;
//	int64* pf = f.data();
//
//	while (pf[d] == 0)
//		d--;
//	return d;
//}

ostream& operator<<(ostream& cout, Poly& f)
{
	int32 d = deg(f);
	int64* pf = f.data();

	printf("[ ");
	for (int32 i = 0; i < d; i++)
		printf("%lld, ", pf[i]);
	printf("%lld ]", pf[d]);

	return cout;
}

ostream& print_vector(int64* f, int32 len)
{
	len--;

	printf("[ ");
	for (int32 i = 0; i < len; i++)
		printf("%lld, ", f[i]);
	printf("%lld ]", f[len]);

	return std::cout;
}


//############################# 分割线 ###############################


void IncompleteNTT::init(int64 p, int64 w, int32 k, int32 h)
{
	this->p = p;
	this->w = w;
	this->k = k;
	this->h = h;
	this->m = (1L << k);
	this->n = this->h * this->m;

	int64 minv, pinv;
	int64 gcd = exgcd(&minv, &pinv, m, p);
	if (gcd != 1)
		ErrorInfo("%s\n", "gcd(2^k, p) != 1");
	this->factor = minv < 0 ? minv + p : minv;
	this->factor_pre = PreCompute(this->factor, p);
	this->one_pre = PreComputeMod(p);

	this->W.resize(this->m + 1); //重复存储：w^{0} = w^{m}
	this->W_pre.resize(this->m + 1);

	int64 *pW = W.data();
	int64 *pW_pre = W_pre.data();

	int64 wi = 1;	//单位根的幂次
	for (int32 i = 0; i <= this->m; i++)
	{
		pW[i] = wi;
		pW_pre[i] = PreCompute(wi, p);
		wi = (wi * w) % p;
	}

	brvlist.resize(k + 1);
	for (int j = 0; j <= k; j++)
	{
		auto& list = brvlist[j];
		int32 len = (1L << j);
		for (int i = 0; i < len; i++)
			list.push_back(brv(i, j));
	}
}




void IncompleteNTT::FwdNTT(NTTRep& F, Poly& f)
{
	int32 d = deg(f);
	if (d >= this->n)
		ErrorInfo("%s\n", "deg(f) > n - 1");

	if (F.size() < this->n)
		F.resize(this->n);

	int64* pF = F.data();
	int64* pf = f.data();

	FwdNTT(pF, pf, f.size());
}

void IncompleteNTT::FwdNTT(int64* pF, int64* pf, int32 flen)
{
	if (pF != pf)
	{
		int32 nn = min(flen, n);
		for (int32 i = 0; i < nn; i++)
			pF[i] = BarrettMod(pf[i], one_pre, p);
	}

	int64* pW = W.data();
	int64* pW_pre = W_pre.data();

	int64 P = 2 * p;
	for (int32 j = 0; j < k; j++)
	{
		int32 blocknum = 1 << j; //对j层分块
		int32 blocksize = this->n >> j; //对j层块大小
		int32 offset = blocksize >> 1; //计算第j+1层时蝴蝶的偏移量(块大小的一半)
		int32* list = brvlist[j].data();

		if (offset >= 4 * h)
			for (int32 i = 0; i < blocknum; i++)
			{
				int64 X, WY;
				int32 s = i * blocksize; //块的起始位置
				int64 ww = pW[list[i] << (k - j - 1)]; //第j层第i个分块使用的单位根，w_{2^{j+1}}^{brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[list[i] << (k - j - 1)];

				int64* pX = pF + (s);
				int64* pY = pF + (s + offset);

				for (int32 k = 0; k < offset; k += 4)
				{
					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;

					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;

					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;

					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;
				}
			}
		else if (offset == 2 * h)
			for (int32 i = 0; i < blocknum; i++)
			{
				int64 X, WY;
				int32 s = i * blocksize; //块的起始位置
				int64 ww = pW[list[i] << (k - j - 1)]; //第j层第i个分块使用的单位根，w_{2^{j+1}}^{brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[list[i] << (k - j - 1)];

				int64* pX = pF + (s);
				int64* pY = pF + (s + offset);

				for (int32 k = 0; k < h; k++)
				{
					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;

					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;
				}
			}
		else
			for (int32 i = 0; i < blocknum; i++)
			{
				int64 X, WY;
				int32 s = i * blocksize; //块的起始位置
				int64 ww = pW[list[i] << (k - j - 1)]; //第j层第i个分块使用的单位根，w_{2^{j+1}}^{brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[list[i] << (k - j - 1)];

				int64* pX = pF + (s);
				int64* pY = pF + (s + offset);

				for (int32 k = 0; k < h; k++)
				{
					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;
				}
			}
	}

	//迭代结束，结果是m=2^k个h长多项式
	//约束范围到[0,2p)
	if (k >= 2)
		for (int32 i = 0; i < n; i += 4)
		{
			*pF -= ((P - *pF) >> 63)&P;
			pF++;

			*pF -= ((P - *pF) >> 63)&P;
			pF++;

			*pF -= ((P - *pF) >> 63)&P;
			pF++;

			*pF -= ((P - *pF) >> 63)&P;
			pF++;
		}
	else
		for (int32 i = 0; i < n; i++)
		{
			*pF -= ((P - *pF) >> 63)&P;
			pF++;
		}
}



void IncompleteNTT::InvNTT(Poly& f, NTTRep& F)
{
	int32 d = deg(F);
	if (d >= this->n)
		ErrorInfo("%s\n", "deg(F) > n - 1");

	if (f.size() < this->n)
		f.resize(this->n);

	int64* pF = F.data();
	int64* pf = f.data();

	InvNTT(pf, pF);
}

void IncompleteNTT::InvNTT(int64* pf, int64* pF)
{
	if (pF != pf)
	{
		for (int32 i = 0; i < n; i++)
			pf[i] = pF[i];
	}

	int64* pW = W.data();
	int64* pW_pre = W_pre.data();

	int64 P = 2 * p;
	for (int32 j = k - 1; j >= 0; j--)
	{
		int32 blocknum = 1 << j; //对j+1层分块
		int32 blocksize = this->n >> j; //第j+1层块大小
		int32 offset = blocksize >> 1; //计算第j层时蝴蝶的偏移量(块大小的一半)
		int32* list = brvlist[j].data();

		if (offset >= 4 * h)
			for (int32 i = 0; i < blocknum; i++)
			{
				int64 XY, T;
				int32 s = i * blocksize; //块的起始位置
				int64 ww = pW[m - (list[i] << (k - j - 1))]; //第j+1层第i个分块使用的单位根，w_{2^{j+1}}^{-brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[m - (list[i] << (k - j - 1))];

				int64* pX = pf + (s);
				int64* pY = pf + (s + offset);

				for (int32 k = 0; k < offset; k += 4)
				{
					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;

					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;

					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;

					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;
				}
			}
		else if (offset == 2 * h)
			for (int32 i = 0; i < blocknum; i++)
			{
				int64 XY, T;
				int32 s = i * blocksize; //块的起始位置
				int64 ww = pW[m - (list[i] << (k - j - 1))]; //第j+1层第i个分块使用的单位根，w_{2^{j+1}}^{-brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[m - (list[i] << (k - j - 1))];

				int64* pX = pf + (s);
				int64* pY = pf + (s + offset);

				for (int32 k = 0; k < h; k++)
				{
					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;

					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;
				}
			}
		else
			for (int32 i = 0; i < blocknum; i++)
			{
				int64 XY, T;
				int32 s = i * blocksize; //块的起始位置
				int64 ww = pW[m - (list[i] << (k - j - 1))]; //第j+1层第i个分块使用的单位根，w_{2^{j+1}}^{-brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[m - (list[i] << (k - j - 1))];

				int64* pX = pf + (s);
				int64* pY = pf + (s + offset);

				for (int32 k = 0; k < h; k++)
				{
					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;
				}
			}
	}

	//迭代结束，再整理一下
	int64 tmp;
	if (k >= 2)
		for (int32 i = 0; i < n; i += 4)
		{
			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;

			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;

			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;

			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;
		}
	else
		for (int32 i = 0; i < n; i++)
		{
			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;
		}
}


void IncompleteNTT::convolution(int64* FG, int64* F, int64* G, int64 w)
{
	for (int32 i = 0; i < h; i++)
	{
		FG[i] = 0;
		int32 s = h + i;
		for (int32 j = 0; j <= i; j++)
			FG[i] = BarrettMod(FG[i] + G[j] * F[i - j], one_pre, p); //冗余，[0,2p)
		for (int32 j = i + 1; j < h; j++)
		{
			int64 tmp = FG[i] + BarrettMod(w * G[j], one_pre, p) * F[s - j];
			FG[i] = BarrettMod(tmp, one_pre, p); //冗余，[0,2p)
		}
	}
}


void IncompleteNTT::NTTMul(NTTRep& FG, NTTRep& F, NTTRep& G)
{
	if (F.size() < n || G.size() < n)
		ErrorInfo("%s\n", "F.size < n or G.size < n");

	if (FG.size() < n)
		FG.resize(n);

	int64* pFG = FG.data();
	int64* pF = F.data();
	int64* pG = G.data();

	NTTMul(pFG, pF, pG);
}

void IncompleteNTT::NTTMul(int64* pFG, int64* pF, int64* pG)
{
	int64* pW = W.data();

	int64 ww;
	int64 r0, r1;
	int64 P = 2 * p;
	int32* list = brvlist[k].data();
	if (h == 1)
		if (k >= 2)
			for (int32 i = 0; i < m; i += 4) //m个h长小多项式
			{
				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;

				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;

				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;

				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;
			}
		else
			for (int32 i = 0; i < m; i++) //m个h长小多项式
			{
				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;
			}
	else if (h == 2)
		if (k >= 2)
			for (int32 i = 0; i < m; i += 4) //m个h长小多项式
			{
				int64 tmp;

				//使用冗余表示，[0,2p)
				ww = pW[list[i]]; //第k层第i个多项式使用的单位根，w_{2^k}^{brv_k(i)}

				tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;

				//使用冗余表示，[0,2p)
				ww = pW[list[i + 1]]; //第k层第i个多项式使用的单位根，w_{2^k}^{brv_k(i)}

				tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;

				//使用冗余表示，[0,2p)
				ww = pW[list[i + 2]]; //第k层第i个多项式使用的单位根，w_{2^k}^{brv_k(i)}

				tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;

				//使用冗余表示，[0,2p)
				ww = pW[list[i + 3]]; //第k层第i个多项式使用的单位根，w_{2^k}^{brv_k(i)}

				tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;
			}
		else
			for (int32 i = 0; i < m; i++) //m个h长小多项式
			{
				//使用冗余表示，[0,2p)
				ww = pW[list[i]]; //第k层第i个多项式使用的单位根，w_{2^k}^{brv_k(i)}

				int64 tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;
			}
	else
		for (int32 i = 0; i < m; i++) //m个h长小多项式
		{
			ww = pW[list[i]]; //第k层第i个多项式使用的单位根，w_{2^k}^{brv_k(i)}
			convolution(pFG, pF, pG, ww);

			pFG += h;
			pF += h;
			pG += h;
		}
}

void IncompleteNTT::Param()
{
	cout << "p = " << p; pn;
	cout << "w = " << w; pn;
	cout << "k = " << k; pn;
	cout << "h = " << h; pn;
	cout << "n = " << n; pn;
	cout << "factor = " << factor; pn;
}


//############################# 分割线 ###############################

void NegacyclicNTT::init(int64 p, int64 w, int32 k, int32 h, bool tag)
{
	if (tag == 0)
		IncompleteNTT::init(p, w, k, 2 * h); //构建父类IncompleteNTT，并初始化
	else
		IncompleteNTT::init(p, w, k + 1, h); //构建父类IncompleteNTT，并初始化

	hh = h;
	hm = m >> 1;
	hn = n >> 1;

	int64 hminv, pinv;
	int64 gcd = exgcd(&hminv, &pinv, hm, p);
	if (gcd != 1)
		ErrorInfo("%s\n", "gcd(2^{k-1}, p) != 1");
	this->factor = hminv < 0 ? hminv + p : hminv;
	this->factor_pre = PreCompute(this->factor, p);
}


void NegacyclicNTT::FwdNTT(NTTRep& F, Poly& f)
{
	int32 d = deg(f);
	if (d >= this->hn)
		ErrorInfo("%s\n", "deg(f) > hn - 1");

	if (F.size() < this->hn)
		F.resize(this->hn);

	int64* pF = F.data();
	int64* pf = f.data();

	FwdNTT(pF, pf, f.size());
}

void NegacyclicNTT::FwdNTT(int64* pF, int64* pf, int32 flen)
{
	if (pF != pf)
	{
		int32 nn = min(flen, hn);
		for (int32 i = 0; i < nn; i++)
			pF[i] = BarrettMod(pf[i], one_pre, p);
	}

	int64* pW = W.data();
	int64* pW_pre = W_pre.data();

	int64 P = 2 * p;
	for (int32 j = 1; j < k; j++) //从j=1层开始迭代
	{
		int32 blocknum = 1 << j; //对j层分块
		int32 blocksize = this->n >> j; //第j层块大小
		int32 offset = blocksize >> 1; //计算第j+1层时蝴蝶的偏移量(块大小的一半)
		int32* list = brvlist[j].data();

		if (offset >= 4 * hh)
			for (int32 i = blocknum >> 1; i < blocknum; i++) //只处理下半部分
			{
				int64 X, WY;
				int32 s = (i - (blocknum >> 1)) * blocksize; //块的起始位置
				int64 ww = pW[list[i] << (k - j - 1)]; //第j层第i个分块使用的单位根，w_{2^{j+1}}^{brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[list[i] << (k - j - 1)];

				int64* pX = pF + (s);
				int64* pY = pF + (s + offset);

				for (int32 k = 0; k < offset; k += 4)
				{
					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;

					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;

					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;

					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;
				}
			}
		else if (offset == 2 * hh)
			for (int32 i = blocknum >> 1; i < blocknum; i++) //只处理下半部分
			{
				int64 X, WY;
				int32 s = (i - (blocknum >> 1)) * blocksize; //块的起始位置
				int64 ww = pW[list[i] << (k - j - 1)]; //第j层第i个分块使用的单位根，w_{2^{j+1}}^{brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[list[i] << (k - j - 1)];

				int64* pX = pF + (s);
				int64* pY = pF + (s + offset);

				for (int32 k = 0; k < hh; k++)
				{
					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;

					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;
				}
			}
		else
			for (int32 i = blocknum >> 1; i < blocknum; i++) //只处理下半部分
			{
				int64 X, WY;
				int32 s = (i - (blocknum >> 1)) * blocksize; //块的起始位置
				int64 ww = pW[list[i] << (k - j - 1)]; //第j层第i个分块使用的单位根，w_{2^{j+1}}^{brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[list[i] << (k - j - 1)];

				int64* pX = pF + (s);
				int64* pY = pF + (s + offset);

				for (int32 k = 0; k < hh; k++)
				{
					//CT蝴蝶(Harvey，输入输出范围[0,4p))
					X = *pX;
					X -= ((P - X) >> 63)&P;
					WY = BarrettMulMod(*pY, ww, ww_pre, p); //Barrett算法
					*pX = X + WY;
					*pY = X + (P - WY);
					pX++;
					pY++;
				}
			}
	}

	//迭代结束，结果是hm=2^{k-1}个h长多项式
	//约束范围到[0,2p)
	if (k >= 2)
		for (int32 i = 0; i < hn; i += 4)
		{
			*pF -= ((P - *pF) >> 63)&P;
			pF++;

			*pF -= ((P - *pF) >> 63)&P;
			pF++;

			*pF -= ((P - *pF) >> 63)&P;
			pF++;

			*pF -= ((P - *pF) >> 63)&P;
			pF++;
		}
	else
		for (int32 i = 0; i < hn; i++)
		{
			*pF -= ((P - *pF) >> 63)&P;
			pF++;
		}
}



void NegacyclicNTT::InvNTT(Poly& f, NTTRep& F)
{
	int32 d = deg(F);
	if (d >= this->hn)
		ErrorInfo("%s\n", "deg(F) > hn - 1");

	if (f.size() < this->hn)
		f.resize(this->hn);

	int64* pF = F.data();
	int64* pf = f.data();

	InvNTT(pf, pF);
}

void NegacyclicNTT::InvNTT(int64* pf, int64* pF)
{
	if (pF != pf)
	{
		for (int32 i = 0; i < hn; i++)
			pf[i] = pF[i];
	}

	int64* pW = W.data();
	int64* pW_pre = W_pre.data();

	int64 P = 2 * p;
	for (int32 j = k - 1; j >= 1; j--) //回到第j=1层
	{
		int32 blocknum = 1 << j; //对j+1层分块
		int32 blocksize = this->n >> j; //第j+1层块大小
		int32 offset = blocksize >> 1; //计算第j层时蝴蝶的偏移量(块大小的一半)
		int32* list = brvlist[j].data();

		if (offset >= 4 * hh)
			for (int32 i = blocknum >> 1; i < blocknum; i++) //只处理下半部分
			{
				int64 XY, T;
				int32 s = (i - (blocknum >> 1)) * blocksize; //块的起始位置
				int64 ww = pW[m - (list[i] << (k - j - 1))]; //第j+1层第i个分块使用的单位根，w_{2^{j+1}}^{-brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[m - (list[i] << (k - j - 1))];

				int64* pX = pf + (s);
				int64* pY = pf + (s + offset);

				for (int32 k = 0; k < offset; k += 4)
				{
					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;

					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;

					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;

					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;
				}
			}
		else if (offset == 2 * hh)
			for (int32 i = blocknum >> 1; i < blocknum; i++) //只处理下半部分
			{
				int64 XY, T;
				int32 s = (i - (blocknum >> 1)) * blocksize; //块的起始位置
				int64 ww = pW[m - (list[i] << (k - j - 1))]; //第j+1层第i个分块使用的单位根，w_{2^{j+1}}^{-brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[m - (list[i] << (k - j - 1))];

				int64* pX = pf + (s);
				int64* pY = pf + (s + offset);

				for (int32 k = 0; k < hh; k++)
				{
					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;

					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;
				}
			}
		else
			for (int32 i = blocknum >> 1; i < blocknum; i++) //只处理下半部分
			{
				int64 XY, T;
				int32 s = (i - (blocknum >> 1)) * blocksize; //块的起始位置
				int64 ww = pW[m - (list[i] << (k - j - 1))]; //第j+1层第i个分块使用的单位根，w_{2^{j+1}}^{-brv_{j+1}(2i)}
				int64 ww_pre = pW_pre[m - (list[i] << (k - j - 1))];

				int64* pX = pf + (s);
				int64* pY = pf + (s + offset);
				for (int32 k = 0; k < hh; k++)
				{
					//GS蝴蝶(Harvey，输入输出范围[0,2p))
					XY = *pX + *pY;
					XY -= ((P - XY) >> 63)&P;
					T = *pX + (P - *pY);
					*pX = XY;
					*pY = BarrettMulMod(T, ww, ww_pre, p); //Barrett算法
					pX++;
					pY++;
				}
			}
	}

	//迭代结束，再整理一下
	int64 tmp;
	if (k >= 2)
		for (int32 i = 0; i < hn; i += 4)
		{
			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;

			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;

			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;

			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;
		}
	else
		for (int32 i = 0; i < hn; i++)
		{
			tmp = BarrettMulMod(*pf, factor, factor_pre, p);
			*pf = tmp - (tmp >= p)*p;
			pf++;
		}
}



void NegacyclicNTT::NTTMul(NTTRep& FG, NTTRep& F, NTTRep& G)
{
	if (F.size() < hn || G.size() < hn)
		ErrorInfo("%s\n", "F.size < n or G.size < n");

	if (FG.size() < hn)
		FG.resize(hn);

	int64* pFG = FG.data();
	int64* pF = F.data();
	int64* pG = G.data();

	NTTMul(pFG, pF, pG);
}

void NegacyclicNTT::NTTMul(int64* pFG, int64* pF, int64* pG)
{
	int64* pW = W.data();

	int64 ww;
	int64 r0, r1;
	int64 P = 2 * p;
	int32* list = brvlist[k].data() + hm;
	if (h == 1)
		if (k >= 2)
			for (int32 i = 0; i < hm; i += 4) //m个h长小多项式
			{
				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;

				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;

				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;

				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;
			}
		else
			for (int32 i = 0; i < hm; i++) //m个h长小多项式
			{
				//使用冗余表示，[0,2p)
				*pFG = BarrettMod(*pF * *pG, one_pre, p);

				pFG++;
				pF++;
				pG++;
			}
	else if (h == 2)
		if (k >= 3)
			for (int32 i = 0; i < hm; i += 4) //hm个h长小多项式
			{
				int64 tmp;

				//使用冗余表示，[0,2p)
				ww = pW[list[i]]; //第k层第2^{k-1}+i个多项式使用的单位根，w_{2^k}^{brv_k(2^{k-1}+i)}

				tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;

				//使用冗余表示，[0,2p)
				ww = pW[list[i + 1]]; //第k层第2^{k-1}+i个多项式使用的单位根，w_{2^k}^{brv_k(2^{k-1}+i)}

				tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;

				//使用冗余表示，[0,2p)
				ww = pW[list[i + 2]]; //第k层第2^{k-1}+i个多项式使用的单位根，w_{2^k}^{brv_k(2^{k-1}+i)}

				tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;

				//使用冗余表示，[0,2p)
				ww = pW[list[i + 3]]; //第k层第2^{k-1}+i个多项式使用的单位根，w_{2^k}^{brv_k(2^{k-1}+i)}

				tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;
			}
		else
			for (int32 i = 0; i < hm; i++) //hm个h长小多项式
			{
				//使用冗余表示，[0,2p)
				ww = pW[list[i]]; //第k层第2^{k-1}+i个多项式使用的单位根，w_{2^k}^{brv_k(2^{k-1}+i)}

				int64 tmp = ww * BarrettMod(pF[1] * pG[1], one_pre, p);
				r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(tmp, one_pre, p);
				//r0 = BarrettMod(pF[0] * pG[0], one_pre, p) + BarrettMod(ww * BarrettMod(pF[1] * pG[1], one_pre, p), one_pre, p);
				r1 = BarrettMod(pF[0] * pG[1], one_pre, p) + BarrettMod(pF[1] * pG[0], one_pre, p);
				pFG[0] = r0 - (((P - r0) >> 63)&P);
				pFG[1] = r1 - (((P - r1) >> 63)&P);

				pFG += 2;
				pF += 2;
				pG += 2;
			}
	else
		for (int32 i = 0; i < hm; i++) //hm个h长小多项式
		{
			ww = pW[list[i]]; //第k层第2^{k-1}+i个多项式使用的单位根，w_{2^k}^{brv_k(2^{k-1}+i)}
			convolution(pFG, pF, pG, ww);

			pFG += h;
			pF += h;
			pG += h;
		}
}
PolyBarrett
#include "PolyBarrett.h"


// h = f / g mod p
void Div(Poly& h, Poly f, Poly g, int64 p)
{
	int df = deg(f);
	int dg = deg(g);

	int64 a = g[dg];
	int64 ainv = inv(a, p);

	int64* pf = f.data();
	int64* pg = g.data();

	for (int i = df; i >= dg; i--)
	{
		int64 b = pf[i];
		int64 c = (b * ainv) % p;
		h[i - dg] = c;
		pf[i] = 0;

		if (c == 0)
			continue;

		int64* ptr = pf + (i - dg);
		for (int j = 0; j < dg; j++)
		{
			ptr[j] -= c * pg[j];
			ptr[j] %= p;
			ptr[j] = ptr[j] < 0 ? ptr[j] + p : ptr[j];
		}
	}
}


void PolyBarrett::init(uint64 p, Poly& g, int32 acc, int32 h)
{
	this->p = p;
	this->dg = deg(g);
	this->acc = acc;
	this->k = dg + acc;
	this->g = g;

	this->N = 2 * dg + acc;	//因为需要计算f*m
	this->H = h;
	this->K = ceil(log2(N / h));
	this->n = pow(2, K)*H;

	uint64 x = (this->n * p*p) >> K;	// P > n*p^2
	do {
		x += 1;
		P = (x << K) + 1;
	} while (!Miller_Rabin(P, 10));	// 2^k | P-1

	int HN = pow(2, K) / 2;
	for (W = 2; W < P; W++)
		if (pow_mod(W, HN, P) == P - 1)	// order(w)=2^K
			break;

	Poly xk(k + 1);
	xk[k] = 1;
	this->m.resize(this->acc + 1);
	Div(this->m, xk, g, p);

	printf("Barrett of Polynomials: acc=%d, dg=%d, p=%lld -> ", acc, dg, p);
	printf("P=%lld, W=%lld, K=%d, H=%d.\n", P, W, K, H);

	this->NTT.init(P, W, K, H);
	this->NTT.FwdNTT(this->M, this->m);
	this->NTT.FwdNTT(this->G, this->g);
}


void PolyBarrett::Mod(int64* h, int64* f, int flen)
{
	Poly t;
	NTTRep T, F(this->n);

	NTT.FwdNTT(F.data(), f, flen);
	NTT.NTTMul(T, F, M);
	NTT.InvNTT(t, T);

	int64* pt = t.data();
	for (int i = 0; i < this->n - this->k; i++)
	{
		pt[i] = pt[i + k] % p;
	}
	for (int i = this->n - this->k; i < this->n; i++)
	{
		pt[i] = 0;
	}

	NTT.FwdNTT(T, t);
	NTT.NTTMul(T, T, G);
	NTT.InvNTT(t, T);

	for (int i = 0; i < this->dg; i++)
	{
		h[i] = (f[i] - pt[i]) % p;
		h[i] += h[i] < 0 ? p : 0;
	}
}
GFX
#include "GFX.h"


//g = f mod (x^n - 1)
void ModXnSubOne(GFX& g, GFX& f, int n)
{
	int fn = f.n;
	int fd = f.d;
	int fp = f.p;

	if (&g != &f)
		for (int i = 0; i < fn; i++)
		{
			int64* p1 = g.rep[i%n].data();
			int64* p2 = f.rep[i].data();
			for (int j = 0; j < fd; j++)
			{
				p1[j] += p2[j];
				p1[j] %= fp;
			}
		}
	else
		for (int i = n; i < fn; i++)
		{
			int64* p1 = g.rep[i%n].data();
			int64* p2 = f.rep[i].data();
			for (int j = 0; j < fd; j++)
			{
				p1[j] += p2[j];
				p1[j] %= fp;
				p2[j] = 0;
			}
		}
}

//g = f mod m(x),p
void Modm(int64* g, int64* f, int flen, Poly& m, uint64 p)
{
	int df;
	for (df = flen - 1; df >= 0; df--)
		if (f[df] != 0)
			break;
	int dm = deg(m);

	if (df < dm)
	{
		for (int i = 0; i <= df; i++)
			g[i] = f[i] % p;
		return;
	}

	int64* pf = new int64[df + 1];
	for (int i = 0; i < df + 1; i++)
		pf[i] = f[i];

	int64 a = m[dm];
	int64 ainv = inv(a, p);
	int64* pm = m.data();

	for (int i = df; i >= dm; i--)
	{
		int64 b = pf[i];
		int64 c = (b * ainv) % p;
		g[i - dm] = c;
		pf[i] = 0;

		if (c == 0)
			continue;

		int64* ptr = pf + (i - dm);
		for (int j = 0; j < dm; j++)
		{
			ptr[j] -= c * pm[j];
			ptr[j] %= p;
			ptr[j] = ptr[j] < 0 ? ptr[j] + p : ptr[j];
		}
	}

	for (int i = 0; i < dm; i++)
		g[i] = pf[i];
	delete pf;
}

ostream& operator<<(ostream& cout, GFX& f)
{
	int deg1 = f.deg();
	printf("[");
	for (int i = 0; i <= deg1; i++)
	{
		int deg2 = deg(f.rep[i]);
		int64* p = f.rep[i].data();

		printf(" [");
		for (int j = 0; j <= deg2; j++)
		{
			printf(" %lld ", p[j]);
		}
		printf("] ");
	}
	printf("]");
	return cout;
}

//c = a+b
void GFX_Calculator::Add(GFX&c, GFX& a, GFX&b)
{
	int d = a.d;
	int n = a.n;
	uint64 p = a.p;

	for (int i = 0; i < n; i++)
	{
		int64* pa = a.rep[i].data();
		int64* pb = b.rep[i].data();
		int64* pc = c.rep[i].data();

		for (int j = 0; j < d; j++)
		{
			pc[j] = pa[j] + pb[j];
			pc[j] -= ((p - pc[j]) >> 63)&p; //if c>p then c = c-p
		}
	}
}

//c = a-b
void GFX_Calculator::Sub(GFX&c, GFX& a, GFX&b)
{
	int d = a.d;
	int n = a.n;
	uint64 p = a.p;

	for (int i = 0; i < n; i++)
	{
		int64* pa = a.rep[i].data();
		int64* pb = b.rep[i].data();
		int64* pc = c.rep[i].data();

		for (int j = 0; j < d; j++)
		{
			pc[j] = pa[j] - pb[j];
			pc[j] += (pc[j] >> 63)&p; //if c<0 then c = c+p
		}
	}
}


//寻找合适的Z_P，以及2^k次单位根w
void GFX_Calculator::findP(int d, int n, uint64 p, uint64& P, uint64& w, int& k, int h)
{
	int N = 2 * ((2 * d)*n);
	k = ceil(log2(N / h));
	int N2 = pow(2, k);	// N2 * h = h*2^k >= N

	uint64 x = (N2*h * p*p) >> k;	// P > n*p^2
	do {
		x += 1;
		P = (x << k) + 1;
	} while (!Miller_Rabin(P, 10));	// 2^k | P-1

	int HN2 = N2 >> 1;
	for (w = 2; w < P; w++)
		if (pow_mod(w, HN2, P) == P - 1)	// order(w)=N2
			break;
}


//c = a*b，GF(p^d)=Z_p[x]/(m(x))
void GFX_Calculator::Mul(GFX&c, GFX& a, GFX&b)
{
	int d = a.d;
	int n = a.n;
	uint64 p = a.p;
	int BlockSize = 2 * d;
	int BlockNum = 2 * n;
	int Len = BlockSize * BlockNum;

	Poly A(Len);
	Poly B(Len);
	Poly C(Len);

	for (int i = 0; i < n; i++)
	{
		int64* pa = a.rep[i].data();
		int64* pb = b.rep[i].data();
		int64* pA = A.data() + (2 * d)*i;
		int64* pB = B.data() + (2 * d)*i;

		for (int j = 0; j < d; j++)
		{
			pA[j] = pa[j];
			pB[j] = pb[j];
		}
	}

	NTTRep A2, B2, C2;

	intt.FwdNTT(A2, A);
	intt.FwdNTT(B2, B);
	intt.NTTMul(C2, A2, B2);
	intt.InvNTT(C, C2);	//C=A*B

	int64* pC = C.data();
	for (int i = 0; i < Len; i++)
		pC[i] %= p;

	/*cout << "A = " << A; pn;
	cout << "B = " << B; pn;
	cout << "C = " << C; pn;*/

	c.zero();
	for (int i = 0; i < BlockNum; i++)
	{
		int64* pC = C.data() + BlockSize * i;
		PB.Mod(c.rep[i].data(), pC, BlockSize);
	}
}


//c = a*b，GF(p^d)=Z_p[x]/(m(x))
void GFX_Calculator::Mul2(GFX&c, GFX& a, GFX&b)
{
	int d = a.d;
	int n = a.n;
	uint64 p = a.p;
	int BlockSize = 2 * d;
	int BlockNum = 2 * n;
	int Len = BlockSize * BlockNum;

	Poly A(Len);
	Poly B(Len);
	Poly C(Len);

	for (int i = 0; i < n; i++)
	{
		int64* pa = a.rep[i].data();
		int64* pb = b.rep[i].data();
		int64* pA = A.data() + (2 * d)*i;
		int64* pB = B.data() + (2 * d)*i;

		for (int j = 0; j < d; j++)
		{
			pA[j] = pa[j];
			pB[j] = pb[j];
		}
	}

	NTTRep A2, B2, C2;

	intt.FwdNTT(A2, A);
	intt.FwdNTT(B2, B);
	intt.NTTMul(C2, A2, B2);
	intt.InvNTT(C, C2);	//C=A*B

	int64* pC = C.data();
	for (int i = 0; i < Len; i++)
		pC[i] %= p;

	/*cout << "A = " << A; pn;
	cout << "B = " << B; pn;
	cout << "C = " << C; pn;*/

	c.zero();
	for (int i = 0; i < BlockNum; i++)
	{
		int64* pC = C.data() + BlockSize * i;
		Modm(c.rep[i].data(), pC, BlockSize, m, p);
	}
}
结果
GFX Mul
m = [ 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 ]

GFX Mul: d=149,n=53,p=2 -> P=147457, W=22, K=14, H=2.
Barrett of Polynomials: acc=149, dg=149, p=2 -> P=3329, W=17, K=8, H=2.
a = [ [ 1 ]  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  [ 1 ] ]

b = [ [ 1  0  1 ]  [ 1 ] ]

10 rounds

使用Barrett：a*b = [ [ 1  0  1 ]  [ 1 ]  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  [ 1  0  1 ]  [ 1 ] ]
0.155162 s

使用长除法：a*b = [ [ 1  0  1 ]  [ 1 ]  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  [ 1  0  1 ]  [ 1 ] ]
0.291891 s

用长除法计算：a*b mod (x^n - 1) = [ [ 0  0  1 ]  [ 1 ]  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  []  [ 1  0  1 ] ]
0.0015902 s
山登绝顶我为峰 3(^v^)3
关注
  • 2
    点赞
  • 踩
  • 6
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
专栏目录
大整数乘法器的FPGA设计与实现
03-15
大整数乘法是公钥加密中最为核心的计算环节，实现运算快速的大数乘法单元是RSA, ElGamal，全同态等密码体制中急需解决的问题之一。针对全同态加密(FHE)应用需求，该文提出一种基于Schönhage-Strassen算法(SSA)的768 kbit大整数乘法器硬件架构。采用并行架构实现了其关键模块64K点有限域快速数论变换(NTT)的运算，并主要采用加法和移位操作以保证并行处理的最大化，有效提高了处理速度。该大整数乘法器在.Stratix-V FPGA上进行了硬件验证，通过与CPU上使用数论库(NTL)和GMP库实现的大整数乘法运算结果对比，验证了该文设计方法的正确性和有效性。实验结果表明，该方法实现的大整数乘法器运算时间比CPU平台上的运算大约有8倍的加速。
有限域GF(2^128)上的乘法，C代码实现
qq_1328512592的博客
05-10 5246
参考：https://blog.csdn.net/bupt073114/article/details/27382533 其中的是GF(2^8)上的乘法，程序、原理讲解很明白了。 一、GF（2^8）上的乘法（直接摘抄自以上网页） //G(28)上的乘法 unsigned char XTIME(unsigned char x) { return ((x << 1) ^ ((x...
有限域GF(2^128)GHASH校验值运算模块——Verilog实现
yinghancidian的博客
09-28 1336
（1）题目详见：有限域GF(2^128)GHASH校验值运算模块说明 （2）实现原理参考的是：有限域GF(2^8)内乘法代码实现以及原理 （3）python验证：有限域GF(2^128)GHASH校验值运算模块——Python验证有限域乘法正确性 （4）仿真结果：有限域GF(2^128)GHASH校验值运算模块——Verilog仿真 这里没有全部按照题目的要求实现算法，而是参照上面超链接（2）中GF(2^8)的实现方式，也就是把乘数B依次左移，然后根据A的每一位的值进行异或。 并行方面，本来想按照题目时序
有限域GF(2^128)GHASH校验值运算模块——Verilog仿真
yinghancidian的博客
09-28 728
一共两个部分（1）有限域乘法模块cal_vz()的仿真（2）顶层模块question_2()的仿真，两个模块的代码都在下面的超链接里，这里就不贴了 Verilog的代码：有限域GF(2^128)GHASH校验值运算模块——Verilog实现 (1)有限域乘法模块cal_vz()的仿真: module cal_vz_sim; // ghash_din ghash_ivin 异或得到A reg[127:0] ghash_din = 128'h5C7A9427_3D9A2AB4_0FD795BD
galois:有限域中的算术和多项式运算
05-17
加洛瓦 有限域中的算术和多项式运算。 安装 $ npm install @guildofweavers/galois --save 例子 import * as galois from '@guildofweavers/galois' ; // create a prime field with a large modulus const field = galois . createPrimeField ( 2n ** 256n - 351n * 2n ** 32n + 1n ) ; const a = field . rand ( ) ; // generate a random field element const b = field . rand ( ) ; // generate another random element const c = field .
在有限域GF(2^n)下求多项式乘法 Python代码实现
weixin_47030183的博客
06-09 4998
#在有限域GF(2^n)下求多项式乘法 Python代码实现 一、了解运算规则 二、例题展示 三、直接上代码（代码有详细备注，不做一一解释） def yxydxscf(a,b,c): # 不可约多项式系数模二运算 c = str(c)[1:] # 不够8位左补零 g = str(a).rjust(8, '0') h = str(b).rjust(8, '0') k = str(c).rjust(8, '0') # 二进制转十进制 d =
Python在GF(2⁸)有限域上求解多项式的乘法逆元——基于扩展欧几里得算法
码疯窝
11-08 1万+
文章目录一、前言二、数学基础1、GF(2⁸)有限域内的多项式2、不可约多项式3、多项式模运算3、乘法逆元三、算法步骤1、扩展欧几里得算法2、多项式除法3、多项式乘法四、代码实现1、多项式除法2、多项式乘法3、欧几里得算法4、配上一个主程序五、心得 一、前言 在写这篇博客前，笔者刚赶完密码学与信息安全的实验报告，用了一天把这个问题涉及的算法仔仔细细啃了一遍，忍不住想写一篇博客。为什么呢？因为笔者昨天...
有限域上乘法
linqing360的博客
04-11 6798
转载自 : http://blog.sina.com.cn/s/blog_4985a6a301000837.html对于有限域GF(256),可以先计算出其乘法表。    在GF(256)中，加法就是异或运算，任意一个元素都可以表示成GF(2)    上的一个最多7次的多项式，    所以    0=000   就是0    1=001   就是1    2=0010就是x+0=x    3=00...
有限域GF(2^8)内乘法代码实现以及原理
热门推荐
bupt073114的专栏
05-28 4万+
在密码学中经常用到有限域的乘法，一般在AES中用到的是GF(2^8)有限域内乘法。什么是有限域呢？有限域通俗的讲就是函数的运算结果全都包含在一个域中，不同于实数域，有限域有一个最大值，所有超过这个最大值的数都会经过一定的方法使他回到这个域中，在密码学中应用很广泛，2^8意味着这个域的最大值是256. 以下代码是GF(2^8)有限域内乘法的C代码实现：
mul.rar_有限域 乘法
09-20
实现有限域中乘法，输入二个普通二级制数，输出在本原多项式的乘法结果
论文研究-加速有限域GF(2m)上Montgomery乘法实现的指令集扩展研究.pdf
07-22
通过扩展多项式乘法指令MULGF2和多项式乘加指令MAGF2来加速Montgomery算法的软件实现。性能分析显示，指令集扩展能够显著提高Montgomery算法的执行效率，特别是同时扩展多项式乘法及乘加指令时效果更佳，且当处理器...
有限域切比雪夫多项式的改进算法 (2011年)
04-25
对计算有限域上切比雪夫多项式的特征多项式算法(CPA)进行了改进,以提高算法的执行速度。首先用蒙哥马利模乘代替原有算法中的普通模乘运算,从而降低单次模乘运算的平均运行时间;其次对蒙哥马利模平方运算的算法流程...
xor_mul.rar_AES_AES VHDL_polynomial_有限域_有限域乘法器
09-14
使用列表法，VHDL语言实现的基于多项式基的有限域乘法器，用于AES算法等对有限域乘法有要求的算法
论文研究 - 在Galois场上搜索具有多项式十进制等效项的Monic不可约多项式
05-30
该方法与在Galois字段GF（pq）上的两个多项式的多项式乘法非常相似，但是执行上有所不同。 多项式的十进制等效项已用于标识基本多项式（BP），EP，IP和可归约多项式（RP）。 单声道RP已通过此方法确定，并已被抵消...
[原]有限域的多项式乘法逆元求解
weixin_30580341的博客
05-15 2054
题目： 求解算法，扩展的欧几里得算法 /* @author tilltheendwjx @blog http://blog.csdn.net/wjh200821或者http://www.cnblogs.com/tilltheendwjx/ */ #include<iostream> using namespace std; int indexofmax1(i...
密码学实验_3_GF有限域上的多项式乘法(python 实现)
dc12499574的博客
06-16 1027
实验
密码学--多项式模乘（有限域GF(2^8)内计算）
最新发布
qq_31804097的博客
10-22 1500
这是期末考试必考榜单的第五项。
gf(2 4)有限域的乘法c语言实现,有限域GF(2^n)的C语言实现浅析
weixin_36349685的博客
05-17 791
由于项目的需要，在网上扒了半天，没有找到域GF(2^n)的C语言实现的系统的介绍。本文试图解释偶特征有限域的实现，让读者不必像我一样浪费太多时间在搜索中。本文以GF(2^8)为例。转载请注明出处，谢谢！甲、有限域的加法实现简单的异或运算即可：unsigned char add(unsigned char a, unsigned char b) {return a ^ b;}乙、有限域的减法实现与加...
扩展欧几里得算法(Bezout恒等式)求有限域上多项式的乘法逆元
超级无敌小胖头
05-22 1万+
前言 最近在复习现代密码理论中的AES，AES中的字节变换的核心操作就是求GF(28)GF(2^8)GF(28)上的多项式逆元，这个问题困扰了我一段时间，今天终于得到解决，其实计算方式和数论中求两个数的Bezout算法是一样的，这里感谢数论老师教给我们的用矩阵行初等变换的方法求Bezout，进而求逆元。 Bezout恒等式 设a,b∈Za,b\in\mathcal{Z}a,b∈Z，则a,ba,ba...
一元多项式加法运算的C语言实现
10-18
好的，以下是一元多项式加法运算的C语言实现： ```c #include <stdio.h> #include <stdlib.h> typedef struct node { int coef; // 系数 int expn; // 指数 struct node *next; } Node, *Polynomial; void CreatePolyn(Polynomial *p, int m) { *p = (Polynomial)malloc(sizeof(Node)); (*p)->next = NULL; Node *tail = *p; for (int i = 0; i < m; i++) { Node *node = (Node*)malloc(sizeof(Node)); scanf("%d %d", &node->coef, &node->expn); tail->next = node; tail = node; } tail->next = NULL; } void PrintPolyn(Polynomial p) { Node *node = p->next; while (node) { printf("%d %d ", node->coef, node->expn); node = node->next; } printf("\n"); } void AddPolyn(Polynomial *pa, Polynomial *pb) { Node *a = (*pa)->next, *b = (*pb)->next, *pre = *pa; while (a && b) { if (a->expn == b->expn) { a->coef += b->coef; if (a->coef == 0) { pre->next = a->next; free(a); a = pre->next; } else { pre = a; a = a->next; } b = b->next; } else if (a->expn < b->expn) { pre->next = b; b = b->next; pre->next->next = a; pre = pre->next; } else { pre = a; a = a->next; } } if (b) { pre->next = b; } free(*pb); } int main() { Polynomial pa, pb; int m, n; scanf("%d", &m); CreatePolyn(&pa, m); scanf("%d", &n); CreatePolyn(&pb, n); AddPolyn(&pa, &pb); PrintPolyn(pa); return 0; } ```

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值