NTT功能与实现

news2025/1/21 6:26:26

NTT的基础功用与拓展功能:

1.evaluate和interpolate

evaluate的本质是选择n个点(假设f(x)的度为n)，计算得到其值，因此根据定义可以直接进行代入计算。为了加快计算的过程选取 w n w_n wn​的幂次(DFT问题即离散傅里叶变换)，使用FFT算法来加快计算过程，将上述方法记作 N T T ( f ) NTT(f) NTT(f)

interpolate的本质是根据n个点值计算得到对应的系数，据此可以列出方程直接求解或者利用矩阵进行求解（根据插值多项式的唯一性，解唯一）。为了加快计算的过程，当点值中的点都为 w n w_n wn​的幂次时，可以使用 F F T − 1 FFT^{-1} FFT−1来进行计算，将上述方法记作 N T T − 1 ( f ) NTT^{-1}(f) NTT−1(f)

注:关于 w n w_n wn​的选择
根据 w n w_n wn​的定义，要求 w n 0 w_n^{0} wn0​, w n 1 w_n^{1} wn1​… w n n − 1 w_n^{n-1} wnn−1​互不相同且 w n n = 1 w_n^{n}=1 wnn​=1，当f(x)的数值需要mod N时，则需要找到 w n n w_n^{n} wnn​同余1模上N。

在这里插入图片描述

FFT算法实现

P3803 【多项式乘法FFT】

#include<iostream>
#include<cmath>
using namespace std;
#define NMAX 10000007
int n, m;
int rev[NMAX];
struct complex {
	//复数类
	double x, y;//x+y*i的格式
	complex(double xx = 0, double yy = 0) {
		x = xx;
		y = yy;
	}
	complex operator + (const complex b) {
		return complex(x + b.x, y + b.y);
	}
	complex operator - (const complex b) {
		return complex(x - b.x, y - b.y);
	}
	complex operator * (const complex b) {
		return complex(x*b.x-y*b.y, x*b.y + y*b.x);
	}
};
struct complex f[NMAX], g[NMAX];//需要在复数域上进行计算
struct complex Wn(int n,int type) {
	//n代表的是等分的分数,type为1代表返回Wn，type为-1代表返回Wn^(-1)
	double Pi = acos(-1.0);
	return complex(cos(2 * Pi / n), type * sin(2 * Pi / n));

};
void test_for_complex() {
	
	complex a = complex(1, 1);
	complex b = complex(2, 2);
	printf("%f+i*%f\n", a.x, a.y);
	printf("%f+i*%f\n", b.x, b.y);
	complex c = a + b;
	printf("%f+i*%f\n", c.x, c.y);
	c = a - b;
	printf("%f+i*%f\n", c.x, c.y);
	c = a * b;
	printf("%f+i*%f\n", c.x, c.y);
}
void FFT(complex *a,int deg,int deg_len,int type) {
	//1.进行比特反转
	for (int i = 0; i < deg; i++)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (deg_len - 1));
	//for (int i = 0; i < deg; i++) cout << rev[i] << " ";
	for (int i = 0; i < deg; i++) {
		if(i<rev[i]) swap(a[i],a[rev[i]]);//注意这里只能交换一次
	}


	//2.进行迭代计算
	for (int m = 2; m <= deg; m <<= 1) {//m代表的是合并后的个数
		//2.1获取原根
		struct complex wn = Wn(m, type);
		for (int k = 0; k < deg; k += m) {//k代表的是待处理组(a[k...k+m-1])的第一个位置
			//2.2对于每一组a[k...k+m/2-1]+a[k+m/2...k+m-1]=a[k...k+m-1]
			complex w = complex(1, 0);
			for (int j = 0; j < m / 2; j++) {//k+j指向a[k...k+m/2-1]
				complex t = w * a[k + j + m / 2];
				complex u = a[k + j];
				a[k + j] = t + u;
				a[k + j + m / 2] = u - t;
				w = wn * w;
			}
		}
	}

}
int main() {
	
	cin >> n >> m;
	for (int i = 0; i < n+1; i++) cin >> f[i].x;
	for (int j = 0; j < m+1; j++) cin >> g[j].x;
	
	//确定等分的分数（由于需要进行加速，所以分数应为2^n的形式）
	int num = 1, len = 0;
	while (num < (n + m+1)) {
		num <<=  1;
		len++;
	}
	FFT(f, num, len, 1);
	FFT(g, num, len, 1);
	for (int i = 0; i < num; i++) {
		f[i] = f[i] * g[i];
	}
	FFT(f, num,len, -1);
	for (int i = 0; i < n + m + 1; i++) cout << int(f[i].x/num + 0.5) << " ";
}

多项式卷积计算

#include<iostream>
#include<cmath>
#include<string>
#include<string.h>
using namespace std;
#define NMAX 102
int n, m;
int rev[NMAX];
struct complex {
	//复数类
	double x, y;//x+y*i的格式
	complex(double xx = 0, double yy = 0) {
		x = xx;
		y = yy;
	}
	complex operator + (const complex b) {
		return complex(x + b.x, y + b.y);
	}
	complex operator - (const complex b) {
		return complex(x - b.x, y - b.y);
	}
	complex operator * (const complex b) {
		return complex(x*b.x-y*b.y, x*b.y + y*b.x);
	}
};
//struct complex f[NMAX], g[NMAX];//需要在复数域上进行计算
struct complex Wn(int n,int type) {
	//n代表的是等分的分数,type为1代表返回Wn，type为-1代表返回Wn^(-1)
	double Pi = acos(-1.0);
	return complex(cos(2 * Pi / n), type * sin(2 * Pi / n));

};
void test_for_complex() {
	
	complex a = complex(1, 1);
	complex b = complex(2, 2);
	printf("%f+i*%f\n", a.x, a.y);
	printf("%f+i*%f\n", b.x, b.y);
	complex c = a + b;
	printf("%f+i*%f\n", c.x, c.y);
	c = a - b;
	printf("%f+i*%f\n", c.x, c.y);
	c = a * b;
	printf("%f+i*%f\n", c.x, c.y);
}
void FFT(complex *a,int deg,int deg_len,int type) {
	//1.进行比特反转
	for (int i = 0; i < deg; i++)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (deg_len - 1));
	//for (int i = 0; i < deg; i++) cout << rev[i] << " ";
	for (int i = 0; i < deg; i++) {
		if(i<rev[i]) swap(a[i],a[rev[i]]);//注意这里只能交换一次
	}


	//2.进行迭代计算
	for (int m = 2; m <= deg; m <<= 1) {//m代表的是合并后的个数
		//2.1获取原根
		struct complex wn = Wn(m, type);
		for (int k = 0; k < deg; k += m) {//k代表的是待处理组(a[k...k+m-1])的第一个位置
			//2.2对于每一组a[k...k+m/2-1]+a[k+m/2...k+m-1]=a[k...k+m-1]
			complex w = complex(1, 0);
			for (int j = 0; j < m / 2; j++) {//k+j指向a[k...k+m/2-1]
				complex t = w * a[k + j + m / 2];
				complex u = a[k + j];
				a[k + j] = t + u;
				a[k + j + m / 2] = u - t;
				w = wn * w;
			}
		}
	}

}

void FFT_new(complex* a, int deg, int deg_len, int type) {
	//1.进行比特反转
	for (int i = 0; i < deg; i++)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (deg_len - 1));
	//for (int i = 0; i < deg; i++) cout << rev[i] << " ";
	for (int i = 0; i < deg; i++) {
		if (i < rev[i]) swap(a[i], a[rev[i]]);//注意这里只能交换一次
	}

	//迭代计算
	for (int m = 2; m <= deg; m <<=1) {
		complex wn = Wn(m, type);//若为逆FTT，则为wn^(-1)
		for (int j = 0; j < deg; j += m) {
			//处理a[j...j+m-1]
			complex w = complex(1, 0);
			for (int k = 0; k < m / 2; k++) {
				//框定左侧为a[j...j+m/2-1]，由j+k游标指向
				complex t = w * a[j + k + m / 2];//旋转因子
				a[j + k + m / 2] = a[j + k] - t;
				a[j + k] = a[j + k] + t;
				w = wn * w;
			}
		}
	}
	
	//逆FTT需要乘上1/n
	if (type == -1) {
		for (int i = 0; i < deg; i++) a[i].x = int(a[i].x / deg + 0.5);
	}
}
int read(string input, complex * f,int start, int end) {
	//从string[start...end]中剥离出多项式系数
	int deg = 0;
	while(input[start] != '(') start++;//过滤掉不必要的空格使得从左括号开始
	double coe = 0, exp = 0;
	for (int i = start+1; i <= end; i++) {
		if (input[i] == ' ') continue;//遇到空格则直接跳过
		else if (input[i] == '+' || input[i] == ')') {
			//cout << coe << " " << exp << endl;
			f[int(exp)].x = coe;
			if (exp > deg) deg = int(exp);
			coe = 0, exp = 0;
		}
		else if (input[i] == 'a') {
			//注意可能存在系数为1的情况
			if (coe == 0) coe = 1;
			i+=2;//跳过^
			while (input[i] != '+' && input[i] != ')') {
				exp = exp * 10 + input[i++] - '0';
			}
			i--;
		}
		else coe = coe * 10 + input[i] - '0';
	}
	return deg + 1;
}
void output(complex* f,int deg) {
	for (int i = deg; i >= 0; i--) {
		if (f[i].x > 0) {
			if (i == 0) cout << int(f[i].x);
			else if (int(f[i].x) != 1)cout << int(f[i].x) << "a^" << i;
			else cout << "a^" << i;
			if (i == 0)cout << endl;
			else cout << "+";
		}
	}
}
int main() {
	string input;
	while (getline(cin, input)) {
		struct complex f[NMAX], g[NMAX];//需要在复数域上进行计算
		int pos =input.find('*');//根据题意，多项式的乘法仅仅只含两项
		int len = strlen(input.c_str());
		if (!(pos < len && pos >= 0)) {
			cout << input << endl;//不含*，则直接输出
			continue;
		}
		int deg_f = read(input, f,0,pos-1);
		int deg_g = read(input, g,pos+1,len-1);
		int deg = 1, deg_len = 0;
		while (deg < (deg_f + deg_g)) deg <<= 1, deg_len++;
		FFT_new(f, deg, deg_len, 1);
		FFT_new(g, deg, deg_len, 1);
		for (int i = 0; i < deg; i++) f[i] = f[i] * g[i];
		FFT_new(f, deg, deg_len, -1);
		output(f,deg);
	}
}

NTT算法实现

void NTT(int* a, int n, int x) {
	//参数设置:a代表待处理的数组，n为度,x代表的是是否是逆NTT
	int len = 0, cn = n;
	while (cn) {
		len++;
		cn >>= 1;
	}
	len--;
	for (RI i = 1; i < n; ++i) {
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
		//cout <<" " << i << " "<< rev[i] << endl;
	}
	//首先进行比特反转拷贝
	for (RI i = 0; i < n; ++i) if (i < rev[i]) swap(a[i], a[rev[i]]);

	for (RI i = 1; i < n; i <<= 1) {//对于每一层
		RI gn = ksm(G, (mod - 1) / (i << 1)); //代表的是旋转因子

		for (RI j = 0; j < n; j += (i << 1)) { //以j到j+(i<<1)为一组，计算j开始的位置
			RI t1, t2, g = 1;
			for (RI k = 0; k < i; ++k, g = 1LL * g * gn % mod) {
				t1 = a[j + k], t2 = 1LL * g * a[j + k + i] % mod; //数组a是如何处理得到的
				a[j + k] = (t1 + t2) % mod, a[j + k + i] = (t1 - t2 + mod) % mod;
			}
		}
	}
	if (x == 1) return;
	int ny = ksm(n, mod - 2);//计算得到n^(-1)
	reverse(a + 1, a + n); //翻转[1...n-1]位，原因在于，求逆代入的是w^0,w^(-1),w^(-2),...w^(-n+1)
	                                                             // w^0,w^(n-1),w(n-2),...w^1
	                                            //而此次计算代入的是w^0,w^1,w^2,...w^(n-1)，因此进行反转即可
	for (RI i = 0; i < n; ++i) a[i] = 1LL * a[i] * ny % mod;
}

2.计算多项式的乘法

问题概述：

计算 C ( x ) = A ( x ) ∗ B ( x ) C(x)=A(x)*B(x) C(x)=A(x)∗B(x)

解决方法1

直接按照手算的方式，展开计算，示例代码如下所示。假设A和B的度为n，则时间复杂度为 O ( n 2 ) O(n^2) O(n2)

a=[1,9]
b=[1,6]
c=[0]*(len(a)+len(b))
mod = 998244353 
for i in range(len(a)):
    for j in range(len(b)):
        c[i+j] =(c[i+j] + a[i]*b[j]) % mod
for i in range(len(c)):
    print(c[i])
#print(1)

解决办法2

  1. 首先估算 C ( x ) C(x) C(x)的度为 d e g c deg_c degc​
  2. 计算得到 d e g deg deg,使得 d e g = 2 i deg=2^i deg=2i，且 d e g > d e g c deg>deg_c deg>degc​
  3. 计算向量 a = N T T ( A , d e g ) a=NTT(A,deg) a=NTT(A,deg)，向量 b = N T T ( B , d e g ) b=NTT(B,deg) b=NTT(B,deg)
  4. 计算向量 c = a ∗ b c=a*b c=a∗b
  5. 向量 C = N T T − 1 ( c , d e g ) C=NTT^{-1}(c,deg) C=NTT−1(c,deg)对应了 C ( x ) C(x) C(x)的各个系数

主体思想：
由于 C ( x ) C(x) C(x)的度为 d e g c deg_c degc​，因此至少需要 d e g c deg_c degc​个点值来推算 C ( x ) C(x) C(x)的系数，因此需要在 A ( x ) A(x) A(x)和 B ( x ) B(x) B(x)上至少取 d e g c deg_c degc​个点值。由于NTT要求 d e g deg deg为 2 n 2^n 2n，因此需要进行第1步和第2步。

#include<iostream>
using namespace std;
int read() {
	int q = 0; char ch = ' ';
	while (ch < '0' || ch>'9') ch = getchar();
	while (ch >= '0' && ch <= '9') q = q * 10 + ch - '0', ch = getchar();
	return q;
}
#define RI register int
const int mod = 998244353, G = 3, N = 2100000;
int n;
int a[N], b[N], c[N], rev[N];
//根据小费马定理,a^(p-1) 同余 1 mod p (p为素数) a^(p-1)=a*a^(p-2)，因此a^(-1) = a^(p-2)，使用快速幂来计算a^(p-2)
int ksm(int x, int y) {
	//快速幂，计算x^y
	int re = 1;
	for (; y; y >>= 1, x = 1LL * x * x % mod) {
		if (y & 1) re = 1LL * re * x % mod;
	}
	return re;
}
void NTT(int* a, int n, int x) {
	//Q:为什么这里仅仅只是考虑了模式的度，而不考虑模式的具体公式
	//参数设置:a代表待检测的数组，n为度,x代表的是是否是逆NTT
	int len = 0, cn = n;
	while (cn) {
		len++;
		cn >>= 1;
	}
	len--;
	for (RI i = 1; i < n; ++i) {
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
		//cout <<" " << i << " "<< rev[i] << endl;
	}
	//首先进行比特反转拷贝
	for (RI i = 0; i < n; ++i) if (i < rev[i]) swap(a[i], a[rev[i]]);

	for (RI i = 1; i < n; i <<= 1) {//对于每一层
		RI gn = ksm(G, (mod - 1) / (i << 1)); //代表的是旋转因子

		for (RI j = 0; j < n; j += (i << 1)) { //以j到j+(i<<1)为一组，计算j开始的位置
			RI t1, t2, g = 1;
			for (RI k = 0; k < i; ++k, g = 1LL * g * gn % mod) {
				t1 = a[j + k], t2 = 1LL * g * a[j + k + i] % mod; //数组a是如何处理得到的
				a[j + k] = (t1 + t2) % mod, a[j + k + i] = (t1 - t2 + mod) % mod;
			}
		}
	}
	if (x == 1) return;
	int ny = ksm(n, mod - 2);//计算得到n^(-1)
	reverse(a + 1, a + n); //翻转[1...n-1]位，原因在于，求逆代入的是w^0,w^(-1),w^(-2),...w^(-n+1)
	                                                             // w^0,w^(n-1),w(n-2),...w^1
	                                            //而此次计算代入的是w^0,w^1,w^2,...w^(n-1)，因此进行反转即可
	for (RI i = 0; i < n; ++i) a[i] = 1LL * a[i] * ny % mod;
}

void test_for_ntt() {
	//使用NTT计算一般的多项式乘法的注意点:
	//1.首先需要估计结果的度，即为所有式子度之和，
	//2.NTT的度要求为2^n，因此需要找到最小的数，满足2^n的形式，同时需要大于估计的结果的度
	//3.NTT计算一般多项式乘的过程为，首先计算对于不同的函数，orz个不同的变量对应的值
	//然后按照计算公式计算得到对应的结果多项式在orz个不同的变量处的取值
	//最后逆NTT变换，得到结果多项式的系数
	int a[20] = { 1, 9};
	NTT(a, 2, 1);
	NTT(a, 2, -1);
	int b[20] = { 1, 6};
	NTT(b, 2, 1);
	NTT(b, 2, -1);
	int c[20];
	int deg = 4; //NTT中只能使用2^n来进行使用
	NTT(a, deg, 1);
	NTT(b, deg, 1);
	for (int i = 0; i < deg; i++) {
		c[i] = 1LL * a[i] * b[i] % mod;
	}
	NTT(c, deg, -1);
	for (int i = 0; i < deg; i++) {
		cout << c[i] << endl;
	}
}
int main()
{
	test_for_ntt();
}

例题

多项式乘法逆
在这里插入图片描述
在这里插入图片描述
注:公式推导过程链接

整体思路为:
为了计算mod x d e g x^{deg} xdeg，先计算 mod x ( d e g + 1 ) / 2 x^{(deg+1)/2} x(deg+1)/2，然后利用公式计算mod x d e g x^deg xdeg。具体来说，若deg为奇数，则计算 m o d x ( d e g + 1 ) / 2 mod x^{(deg+1)/2} modx(deg+1)/2 ,利用公式计算得到 mod x d e g + 1 x^{deg+1} xdeg+1的逆，又因为当a≡b mod x n x^n xn 且n>m时，a≡b mod x m x^m xm，则计算结果等于模 x d e g x^{deg} xdeg的逆；若deg为偶数，则计算mod x d e g / 2 x^{deg/2} xdeg/2，利用公式计算得到mod x d e g x^{deg} xdeg

#include<iostream>
using namespace std;
int read() {
	int q = 0; char ch = ' ';
	while (ch < '0' || ch>'9') ch = getchar();
	while (ch >= '0' && ch <= '9') q = q * 10 + ch - '0', ch = getchar();
	return q;
}
#define RI register int
const int mod = 998244353, G = 3, N = 2100000;
int n;
int a[N], b[N], c[N], rev[N];
//根据小费马定理,a^(p-1) 同余 1 mod p (p为素数) a^(p-1)=a*a^(p-2)，因此a^(-1) = a^(p-2)，使用快速幂来计算a^(p-2)
int ksm(int x, int y) {
	//快速幂，计算x^y
	int re = 1;
	for (; y; y >>= 1, x = 1LL * x * x % mod) {
		if (y & 1) re = 1LL * re * x % mod;
	}
	return re;
}
void NTT(int* a, int n, int x) {
	//Q:为什么这里仅仅只是考虑了模式的度，而不考虑模式的具体公式
	//参数设置:a代表待检测的数组，n为度,x代表的是是否是逆NTT
	int len = 0, cn = n;
	while (cn) {
		len++;
		cn >>= 1;
	}
	len--;
	for (RI i = 1; i < n; ++i) {
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
		//cout <<" " << i << " "<< rev[i] << endl;
	}
	//首先进行比特反转拷贝
	for (RI i = 0; i < n; ++i) if (i < rev[i]) swap(a[i], a[rev[i]]);

	for (RI i = 1; i < n; i <<= 1) {//对于每一层
		RI gn = ksm(G, (mod - 1) / (i << 1)); //代表的是旋转因子

		for (RI j = 0; j < n; j += (i << 1)) { //以j到j+(i<<1)为一组，计算j开始的位置
			RI t1, t2, g = 1;
			for (RI k = 0; k < i; ++k, g = 1LL * g * gn % mod) {
				t1 = a[j + k], t2 = 1LL * g * a[j + k + i] % mod; //数组a是如何处理得到的
				a[j + k] = (t1 + t2) % mod, a[j + k + i] = (t1 - t2 + mod) % mod;
			}
		}
	}
	if (x == 1) return;
	int ny = ksm(n, mod - 2);//计算得到n^(-1)
	reverse(a + 1, a + n); //翻转[1...n-1]位，原因在于，求逆代入的是w^0,w^(-1),w^(-2),...w^(-n+1)
	                                                             // w^0,w^(n-1),w(n-2),...w^1
	                                            //而此次计算代入的是w^0,w^1,w^2,...w^(n-1)，因此进行反转即可
	for (RI i = 0; i < n; ++i) a[i] = 1LL * a[i] * ny % mod;
}
void work(int deg, int* a, int* b) {
	//为了计算mod x^deg，先计算 mod x^((deg+1)/2)，然后利用公式计算mod x^deg 
	//若deg为奇数，则计算mod x^((deg+1)/2) ,利用公式计算得到 mod x^(deg+1)，又因为a==b mod x^n -> a==b mod x^m (n>m) ，则计算结果是可以适用于x^deg的
	//若deg为偶数，则计算mod x^(deg/2)，利用公式计算得到mod x^deg
	//公式的计算过程如下（整个计算过程不涉及mod x^n，而是利用NTT计算得到一般的多项式乘法）
	//使用NTT计算一般的多项式乘法的注意点:
	//1.首先需要估计结果的度，即为所有式子度之和，此处的估计结果为2*deg，即deg<<1
	//2.NTT的度要求为2^n，因此需要找到最小的数，满足2^n的形式，同时需要大于估计的结果的度，此处为orz
	//3.NTT计算一般多项式乘的过程为，首先计算对于不同的函数，orz个不同的变量对应的值，即NTT(c, orz, 1), NTT(b, orz, 1);
	//然后按照计算公式计算得到对应的结果多项式在orz个不同的变量处的取值
	//最后逆NTT变换，得到结果多项式的系数，此处NTT(b, orz, -1);
	if (deg == 1) { b[0] = ksm(a[0], mod - 2); return; }
	work((deg + 1) >> 1, a, b);//这里为什么一定要加上1 a==b mod x^n -> a==b mod x^m (n>m)

	//处理度，取orz为大于2*deg的最小2^n
	RI len = 0, orz = 1;
	while (orz < (deg << 1)) orz <<= 1, ++len;
	cout << deg << " " << orz << endl;
	for (RI i = 1; i < orz; ++i) {
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
		//cout <<" " << i << " "<< rev[i] << endl;
	}
	
	//将数组a复制到数组c，即c同于a mod x^deg
	for (RI i = 0; i < deg; ++i) c[i] = a[i];
	for (RI i = deg; i < orz; ++i) c[i] = 0;

	NTT(c, orz, 1), NTT(b, orz, 1);
	for (RI i = 0; i < orz; ++i)
		b[i] = 1LL * (2 - 1LL * c[i] * b[i] % mod + mod) % mod * b[i] % mod; //算数运算符% * /的优先级是一样的，从左到右计算即可
	NTT(b, orz, -1);
	//输出普通的多项式乘法
	for (RI i = 0; i < orz; ++i) printf("%d,", b[i]);
	cout << endl;

	//进行模x^deg处理
	for (RI i = deg; i < orz; ++i) b[i] = 0; //计算得到b mod x^deg
	for (RI i = 0; i < orz; ++i) printf("%d,", b[i]);
	cout << endl;
}
void test_for_ntt() {
	//使用NTT计算一般的多项式乘法的注意点:
	//1.首先需要估计结果的度，即为所有式子度之和，
	//2.NTT的度要求为2^n，因此需要找到最小的数，满足2^n的形式，同时需要大于估计的结果的度
	//3.NTT计算一般多项式乘的过程为，首先计算对于不同的函数，orz个不同的变量对应的值
	//然后按照计算公式计算得到对应的结果多项式在orz个不同的变量处的取值
	//最后逆NTT变换，得到结果多项式的系数
	int a[20] = { 1, 9};
	NTT(a, 2, 1);
	NTT(a, 2, -1);
	int b[20] = { 1, 6};
	NTT(b, 2, 1);
	NTT(b, 2, -1);
	int c[20];
	int deg = 4; //NTT中只能使用2^n来进行使用
	NTT(a, deg, 1);
	NTT(b, deg, 1);
	for (int i = 0; i < deg; i++) {
		c[i] = 1LL * a[i] * b[i] % mod;
	}
	NTT(c, deg, -1);
	for (int i = 0; i < deg; i++) {
		cout << c[i] << endl;
	}
}
int main()
{
	//test_for_ntt();
	//cout << 3 % 5 * 2 << endl;
	n = read();
	for (RI i = 0; i < n; ++i) a[i] = read();
	work(n, a, b);
	for (RI i = 0; i < n; ++i) printf("%d ", b[i]);
	return 0;
}

3.循环卷积

#include<iostream>
#include<cmath>
#define NMAX 2000
using namespace std;
int rev[NMAX];

int n, m;

struct  complex {
	double x, y;//x+y*i
	complex(double xx = 0, double yy = 0) { x = xx, y = yy; }//构造函数
	complex operator +(const  complex b) {
		return complex(x + b.x, y + b.y);
	}
	complex operator -(const complex b) {
		return complex(x - b.x, y - b.y);
	}
	complex operator *(const complex b) {
		return complex(x * b.x - y * b.y, x * b.y + y * b.x);
	}
}f[NMAX],g[NMAX];
struct complex Wn(int n,int type) {
	double Pi = acos(-1.0);
	return complex(cos(2 * Pi / n), type*sin(2 * Pi / n));
}
void FFT(struct complex* f, int deg, int deg_len,int type) {
	//ftt的逆直接取Wn^(-1)
	//首先进行比特反转
	for (int i = 0; i < deg; i++) {
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (deg_len - 1));
	}
	for (int i = 0; i < deg; i++) {
		if (i < rev[i]) swap(f[i], f[rev[i]]);
	}

	//迭代进行计算
	for (int gap = 2; gap <= deg; gap <<= 1) {
		complex w = Wn(gap,type);
		for (int g_start = 0; g_start < deg; g_start += gap) {
			complex x = complex(1, 0);
			for (int start = g_start; start < g_start + gap / 2; start++) {
				complex u = f[start], v = f[start + gap / 2] * x;
				f[start] = u + v;
				f[start + gap / 2] = u - v;
				x = x * w;
			}
		}
	}

	if (type == -1) {
		for (int i = 0; i < deg; i++) {
			f[i].x = f[i].x/ deg;
		}
	}
}
class NTT {
	int a[NMAX] = { 1,9 }, b[NMAX] = { 1,6 };
	int rev[NMAX] = {0};
	int mod = 998244353;//模数
	int g = 3;//mod简化剩余系上的生成元
	int ksm(int a, int b, int mod) {
		int ret = 1;
		while (b) {
			if (b & 1) ret = ((long long )ret) * a % mod;
			a = ((long long )a) * a % mod;
			b >>= 1;
		}
		return ret;
	}
	void ntt(int* f, int deg, int deg_len, int type) {
		//首先进行比特反转
		for (int i = 0; i < deg; i++) {
			rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (deg_len - 1));
		}
		for (int i = 0; i < deg; i++) {
			if (i < rev[i]) swap(f[i], f[rev[i]]);
		}
		//迭代进行计算
		for (int gap = 2; gap <= deg; gap <<= 1) {
			int w = ksm(g, (mod - 1) / gap, mod);
			for (int g_start = 0; g_start < deg; g_start += gap) {
				int x = 1;
				for (int start = g_start; start < g_start + gap / 2; start++) {
					int u = f[start], v = (1LL * f[start + gap / 2] * x) % mod;
					f[start] = (u + v)% mod;
					f[start + gap / 2] = (u - v + mod) % mod;//做减法要加上mod来避免负数出现
					x = (1LL* x * w) % mod;
				}
			}
		}
		//求逆运算处理
		if (type == -1) {
			for (int i = 1; i < deg / 2; i++) swap(f[i], f[deg - i]);
			int inv_deg = ksm(deg, mod - 2, mod);
			for (int i = 0; i < deg; i++) f[i] = 1LL * f[i] * inv_deg % mod;
		}
		//结果输出展示
		for (int i = 0; i < deg; i++) cout << f[i] << " ";
		cout << endl;
	}
public:
	void test_for_ntt() {
		ntt(a, 4, 2, 1);
		ntt(b, 4, 2, 1);
		for (int i = 0; i < 4; i++) a[i] = (1LL *a[i] * b[i]) % mod;
		ntt(a, 4, 2, -1);
		for (int i = 0; i < 4; i++) cout << a[i] << " ";
	}
};

void test_for_ftt() {
	//计算最高次数分别为n和m的多项式f(x)和g(x)的卷积
	cin >> n >> m;
	for (int i = 0; i < n + 1; i++) cin >> f[i].x;
	for (int j = 0; j < m + 1; j++) cin >> g[j].x;

	//确定等分的分数（由于需要进行加速，所以分数应为2^n的形式）
	int num = 1, len = 0;
	while (num < (n + m + 1)) {
		num <<= 1;
		len++;
	}
	FFT(f, num, len, 1);
	FFT(g, num, len, 1);
	for (int i = 0; i < num; i++) {
		//cout << f[i].x << " " << g[i].x << endl;
		f[i] = f[i] * g[i];
	}
	FFT(f, num, len, -1);
	for (int i = 0; i < n + m + 1; i++) cout << int(f[i].x + 0.5) << " ";
}
int main() {
	NTT t;
	t.test_for_ntt();
}

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/966641.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

CXL.cachemem 简介（背景通道）

CXL.cachemem 简介（背景通道）

&#x1f525;点击查看精选 CXL 系列文章&#x1f525; &#x1f525;点击进入【芯片设计验证】社区&#xff0c;查看更多精彩内容&#x1f525; &#x1f4e2; 声明&#xff1a; &#x1f96d; 作者主页&#xff1a;【MangoPapa的CSDN主页】。⚠️ 本文首发于CSDN&#xff0c…
阅读更多...
无涯教程-JavaScript - TDIST函数

无涯教程-JavaScript - TDIST函数

The TDIST function replaces the T.DIST.2T & T.DIST.RT functions in Excel 2010. 描述 该函数返回学生t分布的百分点(概率)​​,其中数值(x)是t的计算值,将为其计算百分点。 t分布用于小样本数据集的假设检验。使用此函数代替t分布的临界值表。 语法 TDIST(x,deg_fr…
阅读更多...
iOS 设置下载部分文件，如何获取完整文件的大小

iOS 设置下载部分文件，如何获取完整文件的大小

在视频的需求中&#xff0c;遇到这样一个需求&#xff0c;播放一视频的时候&#xff0c;要预下载 后面10条视频&#xff0c;但是只下载后面十条视频的前面1M 实现方法 1 创建请求时设置cacheLength resource [[IdiotResource alloc] init];resource.requestURL task.request…
阅读更多...
屋大人少，凶多吉少

屋大人少，凶多吉少

在这个世界上&#xff0c;包括人在内的万事万物&#xff0c;都是有自己的气场存在的。 那么&#xff0c;人所产生的气场&#xff0c;我们称之为人气。 道法自然&#xff0c;阴阳二象&#xff0c;同样&#xff0c;人的身上也是有阴阳二气&#xff0c; 二气平衡充盈&#xff0c;会…
阅读更多...
【Python】批量下载页面资源

【Python】批量下载页面资源

【背景】 有一些非常不错的资源网站，比如一些MP3资源网站。资源很丰富，但是每一个资源都不大，一个一个下载费时费力，想用Python快速实现可复用的批量下载程序。 【思路】 获得包含资源链接的静态页面，用beautifulsoup分析页面，获得所有MP3资源的实际地址，然后下载。…
阅读更多...
将 Spring Boot 应用程序与 Amazon DocumentDB 集成

将 Spring Boot 应用程序与 Amazon DocumentDB 集成

Amazon DocumentDB&#xff08;与 MongoDB 兼容&#xff09;是一种可扩展、高度持久和完全托管的数据库服务&#xff0c;用于操作任务关键型 MongoDB 工作负载。在 Amazon DocumentDB 上&#xff0c;您可以使用相同的 MongoDB 应用程序代码、驱动程序和工具来运行、管理和扩展工…
阅读更多...
bazel介绍以及其发展历史

bazel介绍以及其发展历史

简介 Bazel Google开源的&#xff0c;是一款与 Make、Maven 和 Gradle 类似的开源构建和测试工具。 它使用人类可读的高级构建语言。Bazel 支持多种语言的项目&#xff0c;可为多个平台构建输出。Bazel支持任意大小的构建目标&#xff0c;并支持跨多个代码库和大量用户的大型代…
阅读更多...
JavaSE基础（1）

JavaSE基础（1）

1 初识Java 知识导图 1.1 Java简介及发展史 1.1.1 Java是什么 Java是一种优秀的程序设计语言&#xff0c;它具有令人赏心悦目的语法和易于理解的语义。 不仅如此&#xff0c;Java还是一个有一系列计算机软件和规范形成的技术体系&#xff0c;这个技术体系提供了完整的用于软…
阅读更多...
系统学习Linux-PXE无人值守装机（附改密）

系统学习Linux-PXE无人值守装机（附改密）

目录 pxe实现系统自动安装pxe工作原理 大致的工作过程如下&#xff1a; PXE的组件&#xff1a; 一、配置vsftpd 二、配置tftp 三、准备pxelinx.0文件、引导文件、内核文件 四、配置dhcp 配置ip 配置dhcp 五、创建default文件 六、新建测试主机用来测试装机效果 七、…
阅读更多...
github指南

github指南

记录一些Github上的宝贵功能 1、github trending/star 1.1查看github的热门趋势 可以选择点击上面热门趋势的链接&#xff0c;从主页点进去的方式如下 在这个页面中&#xff0c;你可以选择language/data range等来搜索到你想要的 1.2 查看github的star排行 比如&#xff1a…
阅读更多...
出货量腰斩，不用中国芯片，美国PC巨头要凉了，苹果成为大赢家

出货量腰斩，不用中国芯片，美国PC巨头要凉了，苹果成为大赢家

市调机构Canalys公布的二季度数据显示美国PC巨头戴尔在中国市场的出货量同比暴跌52%&#xff0c;显示出它在公开宣布舍弃中国芯片之后&#xff0c;中国消费者正纷纷抛弃它&#xff0c;毕竟如今的PC品牌如此之多&#xff0c;完全有更多的替代选择。 Canalys公布的数据显示&#…
阅读更多...
Python Opencv实践 - 拉普拉斯（Laplacian）算子边缘检测

Python Opencv实践 - 拉普拉斯（Laplacian）算子边缘检测

import cv2 as cv import numpy as np import matplotlib.pyplot as pltimg cv.imread("../SampleImages/pomeranian.png", cv.IMREAD_GRAYSCALE) print(img.shape)#拉普拉斯边缘检测 #cv.Laplacian(src, ddepth, dst, ksize, scale, delta, borderType) #src:原图 …
阅读更多...
正中优配：沪指震荡涨0.23%，保险、酿酒等板块走强，半导体板块下挫

正中优配：沪指震荡涨0.23%，保险、酿酒等板块走强，半导体板块下挫

1日早盘&#xff0c;沪指、深成指盘中强势震动上扬&#xff0c;创业板指回落翻绿&#xff0c;科创50指数跌超1%&#xff1b;两市半日成交缺乏5000亿元。 到午间收盘&#xff0c;沪指涨0.23%报3127.19点&#xff0c;深成指涨0.25%&#xff0c;创业板指跌0.23%&#xff0c;科创5…
阅读更多...
Ubuntu22.04安装Mongodb7.0

Ubuntu22.04安装Mongodb7.0

Ubuntu安装Mongodb 1.平台支持2.安装MongoDB社区版2.1导入包管理系统使用的公钥2.2为MongoDB创建列表文件2.3重新加载本地包数据库2.4安装MongoDB包1.安装最新版MongoDB2.安装指定版MongoDB 3.运行MongoDB社区版1.目录2.配置文件3.初始化系统4.启动MongoDB5.验证MongoDB是否成功…
阅读更多...
列表、字典的删除操作

列表、字典的删除操作

1.列表的删除操作&#xff0c;可以使用del 列表[索引]、列表.pop(索引)、列表.remove(元素)、列表.clear() del(如果不指定列表索引&#xff0c;就是删除整个列表&#xff0c;再使用就会出现 name XX is not defined) a [10, 20, 30, 40] print(f"删除前为:{a}") de…
阅读更多...
简易虚拟培训系统-UI控件的应用5

简易虚拟培训系统-UI控件的应用5

目录 Toggle控件简介 示例-使用Toggle组实现主轴速度选择 本篇介绍UI控件Toggle&#xff0c;尝试一个小示例-使用单选框实现速度的选择控制。 Toggle控件简介 1. Toggle的结构如下&#xff1a;最重要的Toggle组件挂在Toggle节点上&#xff0c;下面的Image组件用于显示单选框…
阅读更多...
大数据Flink（七十一）：SQL的时间属性

大数据Flink（七十一）：SQL的时间属性

文章目录 SQL的时间属性 一、Flink三种时间属性简介
阅读更多...
【jsthree.js】全景vr看房进阶版

【jsthree.js】全景vr看房进阶版

three小结&#xff1a; Scene场景 指包含了所有要渲染和呈现的三维对象、光源、相机以及其他相关元素的环境&#xff1b;场景可以被渲染引擎或图形库加载和处理&#xff0c;以生成最终的图像或动画 常见属性&#xff1a; scene.background new THREE.Color(0x000000); // …
阅读更多...
【01背包理论】01背包问题dp[j]（滚动数组） ＜动态规划板子＞

【01背包理论】01背包问题dp[j]（滚动数组） ＜动态规划板子＞

【01背包理论】01背包问题dp[j] 有 n 件物品和一个最多能背重量为 w 的背包。 第 i 件物品的重量是 weight[i]&#xff0c;得到的价值是 value[i] 。 每件物品只有一个&#xff0c;求解将哪些物品装入背包里物品价值总和最大。 题解 动态规划 确定 dp 数组以及下标的含义 滚…
阅读更多...
wxpython： 数字时钟秒表

wxpython： 数字时钟秒表

wxpython 数字时钟秒表&#xff1b;定时器和线程>示例 编写 wx_clock.py 如下 # -*- coding: utf-8 -*- """ 定时器和线程>示例""" import wx import time import threadingclass MyFrame(wx.Frame): def __init__(self):""…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023