博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
多项式相关
阅读量:5070 次
发布时间:2019-06-12

本文共 7325 字,大约阅读时间需要 24 分钟。

省选前准备把多项式搞完。(似乎够折磨人的)

首先FFT和NTT板子请出门右转:
0.加减乘
加减有人用我说吗。
乘就是FFT。
1.多项式求逆
对于一个\(n - 1\)次多项式\(A(x)\),求另一个多项式\(B(x)\),满足\(A(x) * B(x) \equiv 1 \ \ (mod \ \ x ^n)\)。(这里的取模表示只保留小于\(n\)的项的系数)
做法就叫倍增吧。
假设已经求得\(B(x)\),满足\(A(x) * B(x) \equiv 1 \ \ (mod \ \ x ^{\lceil \frac{n}{2} \rceil})\)(我也不知道为啥上取整),要求\(C(x)\)满足\(A(x) * C(x) \equiv 1 \ \ (mod \ \ x ^n)\)
首先可以有\(C(x) \equiv B(x) \ \ (mod \ \ x ^ {\lceil \frac{n}{2} \rceil})\)
然后移项平方得\((B - C) ^ 2 \equiv 0 \ \ (mod \ \ x ^ n)\)。这个觉得挺显然的,因为上面的多项式可以看成最高次项为\(\lceil \frac{n}{2} \rceil - 1\)的一个多项式,只不过系数都是0。然后平方就搞出了一个\(n\)次多项式,系数自然也是0。
接着拆开,两边同乘以A得:\(AB ^ 2 - 2B + C = 0\)
于是\(C = B *(2 - AB)\)
到此为止就做完了。
因为想求\(n\)意义下的就必须先求\(mod \ \ \lceil \frac{n}{2} \rceil\),所以采用递归求解。递归边界就是只有一项,那么\(B(0) = A(0) ^ {mod - 2}\)
其中的乘法用NTT解决。
时间复杂度:\(T(n) = T(\frac{n}{2}) + O(nlogn) = O(nlogn)\)……不会证。

#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
using namespace std;#define enter puts("") #define space putchar(' ')#define Mem(a, x) memset(a, x, sizeof(a))#define In inlinetypedef long long ll;typedef double db;const int INF = 0x3f3f3f3f;const db eps = 1e-8;const ll mod = 998244353;const ll G = 3;const int maxn = 2e6 + 5;inline ll read(){ ll ans = 0; char ch = getchar(), last = ' '; while(!isdigit(ch)) last = ch, ch = getchar(); while(isdigit(ch)) ans = (ans << 1) + (ans << 3) + ch - '0', ch = getchar(); if(last == '-') ans = -ans; return ans;}inline void write(ll x){ if(x < 0) x = -x, putchar('-'); if(x >= 10) write(x / 10); putchar(x % 10 + '0');}int n;int rev[maxn];ll a[maxn], b[maxn], c[maxn];In ll quickpow(ll a, ll b){ ll ret = 1; for(; b; b >>= 1, a = a * a % mod) if(b & 1) ret = ret * a % mod; return ret;}In void ntt(ll* a, int len, int flg){ for(int i = 0; i < len; ++i) if(i < rev[i]) swap(a[i], a[rev[i]]); for(int i = 1; i < len; i <<= 1) { ll gn = quickpow(G, (mod - 1) / (i << 1)); for(int j = 0; j < len; j += (i << 1)) { ll tp1, tp2, g = 1; for(int k = 0; k < i; ++k, g = g * gn % mod) { tp1 = a[j + k], tp2 = g * a[j + k + i] % mod; a[j + k] = (tp1 + tp2) % mod, a[j + k + i] = (tp1 - tp2 + mod) % mod; } } } if(flg == 1) return; int inv = quickpow(len, mod - 2); reverse(a + 1, a + len); for(int i = 0; i < len; ++i) a[i] = a[i] * inv % mod;}In void solve(int deg, ll* a, ll* b) //最后的答案储存在B中,而不是C,这里的C只是充当了一个临时数组 { if(deg == 1) {b[0] = quickpow(a[0], mod - 2); return;} solve((deg + 1) >> 1, a, b); int lim = 0, len = 1; while(len < (deg << 1)) len <<= 1, ++lim; for(int i = 1; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lim - 1)); for(int i = 0; i < deg; ++i) c[i] = a[i]; for(int i = deg; i < len; ++i) c[i] = 0; ntt(c, len, 1), ntt(b, len, 1); for(int i = 0; i < len; ++i) b[i] = (2 * 1LL - c[i] * b[i] % mod + mod) % mod * b[i] % mod; ntt(b, len, -1); for(int i = deg; i < len; ++i) b[i] = 0;}int main(){ n = read(); for(int i = 0; i < n; ++i) a[i] = read(); solve(n, a, b); for(int i = 0; i < n; ++i) write(b[i]), space; enter; return 0;}

2.多项式开根
对于\(A(x)\),找一个多项式\(B(x)\)满足\(B ^ 2 (x) \equiv A(x) \ \ (mod \ \ n)\)
做法和上面很像。
首先有\((B ^ 2 - A) ^ 2 \equiv 0 \ \ (mod \ \ n)\)
然后用初中知识得:\((B ^ 2 + A) ^ 2 \equiv 4AB ^ 2 \ \ (mod \ \ n)\)
移项得\(A = (\frac{B ^ 2 + A}{2B})\)
\(A\)换成\(C ^ 2\)就完事了:\(C = \frac{B ^ 2 + A}{2B}\)
这个除法就用刚学的多项式求逆就好啦。
复杂度还是\(O(nlogn)\)的,因为上面说了,每一层求逆元是\(O(nlogn)\)的。
“递归套递归,复杂度不变”(带劲)

#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
using namespace std;#define enter puts("") #define space putchar(' ')#define Mem(a, x) memset(a, x, sizeof(a))#define In inlinetypedef long long ll;typedef double db;const int INF = 0x3f3f3f3f;const db eps = 1e-8;const int maxn = 2e6 + 5;const ll mod = 998244353;const ll inv2 = 499122177;const ll G = 3;inline ll read(){ ll ans = 0; char ch = getchar(), last = ' '; while(!isdigit(ch)) last = ch, ch = getchar(); while(isdigit(ch)) ans = (ans << 1) + (ans << 3) + ch - '0', ch = getchar(); if(last == '-') ans = -ans; return ans;}inline void write(ll x){ if(x < 0) x = -x, putchar('-'); if(x >= 10) write(x / 10); putchar(x % 10 + '0');}int n, rev[maxn];ll a[maxn], b[maxn], c[maxn], tp[maxn];In ll quickpow(ll a, ll b){ ll ret = 1; for(; b; b >>= 1, a = a * a % mod) if(b & 1) ret = ret * a % mod; return ret;}In void ntt(ll* a, int len, bool flg){ for(int i = 0; i < len; ++i) if(i < rev[i]) swap(a[i], a[rev[i]]); for(int i = 1; i < len; i <<= 1) { ll ng = quickpow(G, (mod - 1) / (i << 1)); for(int j = 0; j < len; j += (i << 1)) { ll g = 1; for(int k = 0; k < i; ++k, g = g * ng % mod) { ll tp1 = a[k + j], tp2 = a[k + j + i] * g % mod; a[k + j] = (tp1 + tp2) % mod; a[k + j + i] = (tp1 - tp2 + mod) % mod; } } } if(flg) return; ll inv = quickpow(len, mod - 2); reverse(a + 1, a + len); for(int i = 0; i < len; ++i) a[i] = a[i] * inv % mod;}In void sol_inv(ll* a, ll* b, int deg){ if(deg == 1) {b[0] = quickpow(a[0], mod - 2); return;} sol_inv(a, b, (deg + 1) >> 1); int len = 1, lim = 0; while(len < (deg << 1)) len <<= 1, ++lim; for(int i = 0; i < deg; ++i) tp[i] = a[i]; for(int i = deg; i < len; ++i) tp[i] = 0; for(int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lim - 1)); ntt(tp, len, 1); ntt(b, len, 1); for(int i = 0; i < len; ++i) b[i] = b[i] * (2 * 1LL - tp[i] * b[i] % mod + mod) % mod; ntt(b, len, 0); for(int i = deg; i < len; ++i) b[i] = 0;}In void sol_sqrt(ll* a, ll* b, int deg){ if(deg == 1) {b[0] = 1; return;} sol_sqrt(a, b, (deg + 1) >> 1); for(int i = 0; i < (deg << 1); ++i) c[i] = 0; sol_inv(b, c, deg); int len = 1, lim = 0; while(len < (deg << 1)) len <<= 1, ++lim; for(int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lim - 1)); ntt(b, len, 1); for(int i = 0; i < len; ++i) b[i] = b[i] * b[i] % mod; ntt(b, len, 0); for(int i = 0; i < deg; ++i) b[i] = (b[i] + a[i]) % mod; //循环不是到len! ntt(b, len, 1); ntt(c, len, 1); for(int i = 0; i < len; ++i) b[i] = b[i] * c[i] % mod * inv2 % mod; ntt(b, len, 0); for(int i = deg; i < len; ++i) b[i] = 0;}int main(){ n = read(); for(int i = 0; i < n; ++i) a[i] = read(); sol_sqrt(a, b, n); for(int i = 0; i < n; ++i) write(b[i]), space; enter; return 0;}

3.多项式求导
(原谅我不会求导,现学的)
首先得知道幂函数求导:\((x ^ a)' = ax ^ {a - 1}\)
所以每一项的导数就是\((ax ^ b)' = abx ^ {b - 1}\)
每一项的导数相加就是多项式的导数(简单不)。

In void get_der(ll* a, ll* b, int n)    //导数英文derivation,特意查的 {    for(int i = 1; i < n; ++i) b[i - 1] = a[i] * i % mod;    b[n - 1] = 0;}

4.多项式积分
(积分就是谁的导数是我——送给像我一样不会积分的人)
知道导数怎么求了,我们就可以反推积分(仔细想想):\(\int ax ^ b = \frac{a}{b + 1} x ^ {b + 1}\)
同理每一项加起来就是多项式积分

In void get_int(ll* a, ll* b, int n)    //积分叫integral啦{    for(int i = 1; i < n; ++i) b[i] = a[i - 1] * inv[i] % mod;    b[0] = 0;}

5.多项式取对数
直接上公式:\(In A = \int \frac{A'}{A}\)
也就是先求导再乘以逆元最后积分一下。

ll tp1[maxn], tp2[maxn];In void get_In(ll* a, ll* b, int n){    get_der(a, tp1, n); sol_inv(n, a, tp2);    int len = 1, lim = 0;    while(len < (n << 1)) len <<= 1, ++lim;    for(int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lim - 1));    ntt(tp1, len, 1); ntt(tp2, len, 1);    for(int i = 0; i < len; ++i) tp1[i] = tp1[i] * tp2[i] % mod;    ntt(tp1, len, 0);    get_int(tp1, b, n);}

转载于:https://www.cnblogs.com/mrclr/p/10367006.html

你可能感兴趣的文章
[kmp]HDU1711 Number Sequence
查看>>
Unity Ioc 学习笔记1
查看>>
并不对劲的多项式求逆
查看>>
c/c++关键字:static
查看>>
怎么通过浏览器下载自己喜欢的音频视频
查看>>
Oulipo HDU 1686 KMP模板
查看>>
Linux学习笔记之Linux启动引导过程
查看>>
tomcat9 gzip
查看>>
几种常见的架构模式
查看>>
对数据库中的表或字段重命名
查看>>
Alpha阶段项目总结
查看>>
正则表达式(进阶篇)
查看>>
你知道的和不知道的sass
查看>>
#ifndef/#define/#endif以及#if defined/#else/#endif使用详解
查看>>
实例分析如何精确C#日期格式到毫秒,DateTime调用
查看>>
小记 百度地图 soso地图 经纬度偏移
查看>>
【C】由printf("%d\t%d\t%d\n",a,a+=(a++),a);引起的思考
查看>>
Eclipse与Android源码中ProGuard工具的使用
查看>>
阿里云服务器ecs配置之安装mysql
查看>>
POJ 2484 A Funny Game【博弈】
查看>>