码迷,mamicode.com
首页 > 其他好文 > 详细

多项式的高级运算

时间:2020-07-27 16:01:58      阅读:76      评论:0      收藏:0      [点我收藏+]

标签:read   函数   span   void   简单   lin   lld   转化   register   

基础:FFT与NTT

多项式求乘法逆元

【模板】多项式乘法逆

技术图片

\(Code:\)

int n;
ll A[N], B[N], C[N], r[N];
ll limi, l;
inline ll quickpow(ll x, ll k)...
inline void ntt(ll *a, int type) {...//此处已经让type = 1的乘inv了
void sol(int deg, ll *a, ll *b) {//b is 逆元
	if (deg == 1) {
		b[0] = quickpow(a[0], P - 2);
		return ;
	}
	sol((deg + 1) >> 1, a, b);
	
	limi = 1, l = 0;
	while (limi <= (deg << 1))	limi <<= 1, ++l;
	for (register int i = 1; i <= limi; ++i)
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
	for (register int i = 0; i < deg; ++i)	C[i] = a[i];//转移到C防变化 
	for (register int i = deg; i < limi; ++i)	C[i] = 0;//多次清空更保险 
	ntt(b, 1); ntt(C, 1);
	for (register int i = 0; i < limi; ++i)//B = 2B‘ - AB‘B‘ = B‘(2 - AB‘)
		b[i] = ((2ll - b[i] * C[i]) % P + P) % P * b[i] % P;
	ntt(b, -1);
	for (register int i = deg; i <= limi; ++i)//多次清空更保险 
		b[i] = 0;
}
int main() {
	read(n);
	for (register int i = 0; i < n; ++i)	read(A[i]);
	sol(n, A, B);
	for (register int i = 0; i < n; ++i)
		printf("%lld ", B[i]);
	return 0;
}

这里给一个简化的板子,供复习:

inline void get_inv(ll *a, ll *b, int deg) {//deg:项数 
	if (deg == 1)	return b[0] = inv(a[0]), void() ; 
    
	get_inv(a, b, (deg + 1) >> 1);
	limi = 1; len = 0;
	while (limi <= (deg << 1)) limi <<= 1, len++;
	for (register int i = 0; i < limi; ++i)
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (len - 1));
        
	for (register int i = 0; i < deg; ++i)	c[i] = a[i], d[i] = b[i];
	for (register int i = deg; i < limi; ++i)	c[i] = d[i] = 0;
    
	ntt(c, 1), ntt(d, 1);
	for (register int i = 0; i < limi; ++i)	c[i] = c[i] * d[i] % P * d[i] % P;
	ntt(c, -1);
	for (register int i = 0; i < deg; ++i)	b[i] = ((2 * b[i] - c[i]) % P + P) % P;
}

多项式开根号

题意

  • 给A(x),其中a0 = 1,求B(x),使得B(x)^2 = A(x) (mod x^n)

思路简析

与多项式求逆相同,由于一平方就mod x^m -> mod x^(2m),我们考虑递归求解。

表达式

同样假设我们已经求出B(x)的一半b(x),那么:

A * b = 1(mod x^m)
A * B = 1(mod x^m)
∴B - b = 0(mod x^m)
两边平方:
B^2 - 2 * B * b + b^2 = 0(mod x^(2*m))
据B^2 = A(mod x^(2*m)):
A - 2 * B * b + b^2 = 0(mod x^(2*m))
于是:
B = (A/b + b) / 2

配合多项式求逆解出B(x)。

递归边界

模板题非常友善,告诉我们 a0 = 1,于是递归边界为

if (deg == 1) {b[0] = 1; return ;}

如果题目没有那么友善,那么我们或许可以多random几个数 我们需要用二次剩余之类的麻烦的东西,或者考虑换一种算法。

浅谈二次剩余

Code:

void get_sqrt(ll *a, ll *b, int deg) {
	if (deg == 1) {b[0] = 1; return ;}
	get_sqrt(a, b, (deg + 1) >> 1);
  //get_len
	limi = 1, len = 0;
	while (limi <= (deg << 1))	limi <<= 1, len++;
	for (register int i = 0; i <= limi; ++i)
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (len - 1));
        
  //copy and multiply
	for (register int i = 0; i <= limi; ++i)	bn[i] = 0;//attention
	get_inv(b, bn, deg);
	for (register int i = 0; i < deg; ++i)	C[i] = a[i];
	for (register int i = deg; i <= limi; ++i)	C[i] = 0;
	ntt(C, 1); ntt(bn, 1);
	for (register int i = 0; i <= limi; ++i)
		C[i] = C[i] * bn[i] % P;
	ntt(C, -1);
	for (register int i = 0; i < deg; ++i)	b[i] = (C[i] + b[i]) * inv2 % P;
	for (register int i = deg; i <= limi; ++i)	b[i] = 0;
}

注意:

  1. 用数组前一定注意清空。我也不知道为什么,反正不清空就会出错。估计是NTT的祸吧。

多项式除法

  • P4512 【模板】多项式除法

  • A(x) * B(x) + C(x) = D(x),给出D(x), A(x),求B(x),C(x)。(类似高精除)

  • m = deg(A) < 100,000, n = deg(D) <= 100,000,m <= n

思路简析

我们发现神奇的事情:把A,B,C,D都翻转过来,(把C加到D后面),等式仍然成立。并且还有一个好处:C转过来后D的0~n-m项都不受C的影响,而B又肯定超不过n-m+1项。因此我们可以借助反转后的A,D数组算出B数组,然后什么都好搞了。

Code:

int main() {
	read(n); read(m); n++; m++;//n,m变成项数
	for (register int i = 0; i < n; ++i)	read(D[i]), Dbp[i] = D[i];
	for (register int i = 0; i < m; ++i)	read(A[i]), Abp[i] = A[i];//backup
	Reverse(A, m);
	Reverse(D, n);
	for (register int i = n - m + 1; i < n; ++i)
		A[i] = D[i] = 0;
	get_inv(A, An, n - m + 1);
	mul(D, An, B, n - m + 1, n - m + 1);
    //D(n-m+1项) * An(n-m+1项) -> B
	Reverse(B, n - m + 1);
	mul(Abp, B, AB, m, n - m + 1);
	for (register int i = 0; i < m - 1; ++i)
		C[i] = (Dbp[i] - AB[i] + P) % P;
    ...
}

注意

  • 在算D*inv(A)时,保险起见,只保留D和A的0~n-m项,且对A,D做备份,算C时用。

  • 什么时候用n-m,什么时候用n-m+1,要分清楚!

  • 此时多项式变量名逐渐增多,注意区分,不要把An写成A!

分治FFT

(其实我觉得对于我这个几乎只写NTT的人来说,叫分治NTT比较好)

【模板】分治 FFT

简单说一下,分治FFT用到了CDQ分治的思想,但不用非得学CDQ分治,毕竟这个思想还是比较好理解的,之前也经常用到。简而言之,这里的CDQ分治思想就是:每次只考虑左半段对右半段的贡献,先递归解决左半段,然后让右半段加上左半段的贡献,再递归解决右半段。这样,一次次贡献的加和就组成了每个位置的值。

将题意稍稍转化,f[i] = g[0 ~ i]与f[0 ~ i]的卷积(多项式乘法)。剩下的推导部分可以通过手玩来推导。

最关键的两条语句:

    for (int i = l;i <= mid; i++) a[i-l] = f[i];
    for (int i = 0;i < len; i++) b[i] = g[i];

\(Code:\) my code

(主要篇幅是NTT,其实重点就在于以上两条语句,NTT只是工具)

多项式求ln

在学习微积分后,我再学ln,感觉舒适了很多。

求ln很简单,两边求个到,用一下多项式求逆,再积分即可。

思维难度低,代码量大。

新增两条背诵语句:

inline void dao(ll *a, ll *b, int n) {
	for (register int i = 1; i < n; ++i)	b[i - 1] = i * a[i] % P; b[n - 1] = 0;
}
inline void ji(ll *a, ll *b, int n) {
	for (register int i = 1; i < n; ++i)	b[i] = a[i - 1] * inv(i) % P; b[0] = 0;
}

其实也不是很难背,考虑 求导/积分 完后的b[i](或b[i - 1])是从哪里转移来的,就很容易理解了。

剩余部分代码:

	get_inv(A, An, n); dao(A, Ad, n);
	
	limi = 1, len = 0;
	while (limi <= (n << 1))	limi <<= 1, len++;
	for (register int i = 0; i <= limi; ++i) 
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (len - 1));
	
	ntt(Ad, 1); ntt(An, 1);
	for (register int i = 0; i <= limi; ++i)	A[i] = Ad[i] * An[i] % P;
	ntt(A, -1);
	ji(A, C, n);
	for (register int i = 0; i < n; ++i)
		printf("%lld ", C[i] % P);

多项式的高级运算

标签:read   函数   span   void   简单   lin   lld   转化   register   

原文地址:https://www.cnblogs.com/JiaZP/p/13384622.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!