多项式半家桶
多项式全家桶是什么?
1e6多点求值?
5e5快速插值?
很不幸,下面的板子里一个都没有qwq
多项式乘法NTT
FFT是什么?我还真不太会qwq
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int P = 998244353;
const int G = 3;
const int G1 = 332748118;
const int MAXN = 3e6 + 7;
int n, m;
ll A[MAXN], B[MAXN];
int L, len, R[MAXN];
inline ll add(ll x, ll y) {
x += y;
if(x >= P)x -= P;
return x;
}
inline ll ksm(ll x, ll y) {
ll ans = 1;
while(y) {
if(y & 1)ans = ans * x % P;
x = x * x % P;
y >>= 1;
}
return ans;
}
inline void init() {
len = 1;
while(len <= n + m)len <<= 1, L++;
for(int i = 0; i < len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
return ;
}
inline void NTT(ll *F, int len, int typ) {
for(int i = 0; i < len; ++i) {
if(i < R[i])swap(F[i], F[R[i]]);
}
for(int mid = 1; mid < len; mid <<= 1) {
ll wn = ksm(typ == 1 ? G : G1, (P - 1) / (mid << 1));
for(int j = 0; j < len; j += (mid << 1)) {
ll w = 1;
for(int k = 0; k < mid; ++k, w = w * wn % P) {
ll x = F[j + k], y = F[j + k + mid] * w % P;
F[j + k] = add(x, y);
F[j + k + mid] = add(x, P - y);
}
}
}
if(typ == -1) {
ll inv = ksm(len, P - 2);
for(int i = 0; i < len; ++i) {
F[i] = F[i] * inv % P;
}
}
}
int main() {
scanf("%d%d", &n, &m);
for(int i = 0; i <= n; ++i) {
scanf("%lld", &A[i]);
}
for(int i = 0; i <= m; ++i) {
scanf("%lld", &B[i]);
}
init();
NTT(A, len, 1);
NTT(B, len, 1);
for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
NTT(A, len, -1);
for(int i = 0; i <= n + m; ++i) {
printf("%lld ", A[i]);
}
return 0;
}
代码分析:
首先是init,注意到我们按顺序合并的序列实际上是二进制翻转后的下标序列,所以对于这个R,我们有其在意义下的翻转二进制位,递推方式也很简单,先把之前的平移下来,然后再把最后一位提到处
其次是NTT过程,其实框架按照单位根反演的式子背诵就好了,因为那个式子连除n个乘上单位根逆元啥的都有....
然后是三层循环,这个很好看,就是我们从长度为1到长度为2n的log层挨个做,然后每层有$2^n/2^g$段,每一段都要和下一段的同相对位置做乘法,而且每一层的某段都要乘上一个wn(段数-1),wn即单位根,相当于线段树的每一层长度之和,就是的
注意循环卷积要开至少两倍空间
【模板】分治 FFT
生成函数做法:
序列f的OGF为F(x)其中f的第1项是1,g的生成函数为G(x),其中g的第0项是0
有
于是可以解得多项式求逆即可
分治FFT做法:
CDQ:Yeah!
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int P = 998244353;
const int G1 = 3;
const int G2 = (P + 1) / 3;
const int MAXN = 4e5 + 7;
int n, L, len;
int F[MAXN], G[MAXN], R[MAXN];
ll A[MAXN], B[MAXN];
inline void init(int x) {
len = 1;
L = 0;
while(len <= 2 * x)len <<= 1, ++L;
for(int i = 1; i <= len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
return ;
}
inline ll ksm(ll x, ll y) {
ll ans = 1;
while(y) {
if(y & 1)ans = ans * x % P;
x = x * x % P;
y >>= 1;
}
return ans;
}
inline void add(int &x, ll y) {
x += y;
if(x >= P)x -= P;
}
inline void NTT(ll *F, int len, int typ) {
for(int i = 1; i < len; ++i)if(i < R[i])swap(F[i], F[R[i]]);
for(int mid = 1; mid < len; mid <<= 1) {
ll wn = ksm(typ == 1 ? G1 : G2, (P - 1) / (mid << 1));
for(int j = 0; j < len; j += (mid << 1)) {
ll w = 1;
for(int k = 0; k < mid; ++k, w = w * wn % P) {
int x = F[j + k], y = F[j + k + mid] * w % P;
F[j + k] = (x + y) % P;
F[j + k + mid] = (x + P - y) % P;
}
}
}
if(typ == -1) {
ll inv = ksm(len, P - 2);
for(int i = 0; i < len; ++i)F[i] = F[i] * inv % P;
}
return ;
}
inline void out1() {
puts("A->:");
for(int i = 0; i <= len; ++i) {
printf("%lld ", A[i]);
}
puts("");
return ;
}
inline void out2() {
puts("B->:");
for(int i = 0; i <= len; ++i) {
printf("%lld ", B[i]);
}
puts("");
return ;
}
inline void cdq(int l, int r) {
if(l == r)return ;
int mid = (l + r) >> 1;
cdq(l, mid);
for(int i = 0; i < r - l + 1; ++i)A[i] = G[i]; //0~r-l+1,注意我们是多项式乘法,所以这个必须是前...项
for(int i = l; i <= mid; ++i)B[i - l] = F[i];
for(int i = mid + 1; i <= r; ++i)B[i - l] = 0;
init(r - l + 1);//传入项数
NTT(A, len, 1);
NTT(B, len, 1);
for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
NTT(A, len, -1);
for(int i = mid + 1; i <= r; ++i)add(F[i], A[i - l]);
for(int i = 0; i < len; ++i)A[i] = B[i] = 0;
cdq(mid + 1, r);
return ;
}
inline void zhqwq() {
for(int i = 0; i < 10; ++i)A[i] = 1;
for(int j = 0; j < 10; ++j)B[j] = 1;//1+x+x^2.....,等价于求前缀和
init(10);
NTT(A, len, 1);
NTT(B, len, 1);
for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
NTT(A, len, -1);
return ;
}
int main() {
scanf("%d", &n);
for(int i = 1; i < n; ++i) {
scanf("%d", &G[i]);
}
init(n);
F[0] = 1;
cdq(0, n - 1);
for(int i = 0; i < n; ++i) {
printf("%d ", F[i]);
}
return 0;
}
代码分析:调试NTT可以使用zhqwq函数,因为只要生成函数学的好,乘上啥得出啥很容易知道,所以NTT就不会错qwq
分治过程我们要注意是拿出左半部分的F,右半部分都是0,而且g数组要用从0~r-l项,因为我们是用前半部分进行卷积!
然后记得清空要把len范围都清空,因为我们是循环卷积,后面不清空会错的
【模板】多项式乘法逆
证明:
如果我们mod的不是x,就把这个mod项减少一下
假设已知
显然的
做差有
于是我们除掉这个F,然后同时平方有
然后我们发现,无论如何左边都会有一项的系数小于所以就会变没了qwq
所以说这个就做完啦sth=0!
那么
同时乘这个,因为
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN = 4e5 + 7;
const int P = 998244353;
const int G1 = 3;
const int G2 = (P + 1) / 3;
int N, len, L, R[MAXN];
ll G[MAXN], A[MAXN], F[MAXN], B[MAXN];
inline void init(int x) {
len = 1;
L = 0;
while(len <= 2 * x)len <<= 1, ++L;
for(int i = 1; i < len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
return ;
}
inline ll ksm(ll x, ll y) {
ll ans = 1;
while(y) {
if(y & 1)ans = ans * x % P;
x = x * x % P;
y >>= 1;
}
return ans;
}
inline void NTT(ll *F, int len, int typ) {
for(int i = 1; i < len; ++i)if(i < R[i])swap(F[i], F[R[i]]);
for(int mid = 1; mid < len; mid <<= 1) {
ll wn = ksm(typ == 1 ? G1 : G2, (P - 1) / (mid << 1));
for(int j = 0; j < len; j += (mid << 1)) {
ll w = 1;
for(int k = 0; k < mid; ++k, w = w * wn % P) {
int x = F[j + k], y = F[j + k + mid] * w % P;
F[j + k] = (x + y) % P;
F[j + k + mid] = (x + P - y) % P;
}
}
}
if(typ == -1) {
ll inv = ksm(len, P - 2);
for(int i = 0; i < len; ++i)F[i] = F[i] * inv % P;
}
return ;
}
inline void out1() {
printf("A-:");
for(int i = 0; i < len; ++i)printf("%d ", A[i]);
puts("\n");
}
inline void getinv(ll *F, int n) {
if(n == 1) {
G[0] = ksm(F[0], P - 2);
return ;
}//只有一项qwq
int m = (n + 1) / 2;
getinv(F, m);
init(n);
for(int i = 0; i < n; ++i)A[i] = G[i], B[i] = F[i];
for(int i = n; i < len; ++i)A[i] = B[i] = 0;
NTT(A, len, 1);
NTT(B, len, 1);
for(int i = 0; i < len; ++i)A[i] = ( 2 * A[i] % P - A[i] * B[i] % P * A[i] % P + P) % P;
NTT(A, len, -1);
for(int i = 0; i < n; ++i)G[i] = A[i];//还原G
return ;
}
inline void zhqwq() {
for(int i = 0; i < 10; ++i)A[i] = 1;
for(int j = 0; j < 10; ++j)B[j] = 1;//1+x+x^2.....,等价于求前缀和
init(10);
NTT(A, len, 1);
NTT(B, len, 1);
for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
NTT(A, len, -1);
return ;
}
int main() {
scanf("%d", &N);
for(int i = 0; i < N; ++i)scanf("%lld", &F[i]);
getinv(F, N);//递归这个mod x^n
for(int i = 0; i < N; ++i) {
printf("%lld ", G[i]);
}
return 0;
}
代码分析:写的时候NTT的w*wn没有取模...服了
然后小心这个getinv我们递归下去是取一半上取整递归哈qwq以及最后到底层是n=1,我们直接求逆即可
【模板】多项式对数函数(多项式 ln)
说到ln,很快啊!,就有chz上课不会初中数学的丢人历史
证明:
求导,链式法则
所以我们直接求出F'和F逆,然后乘一下,然后再积分积回去就好了
这意味我们没有这个常数项qwq,而你会发现这个常数项是1?
取ln后的多项式常数项是0....
所以其他有理数不行也很好证明了,因为e是超越数,所以除了c=0是1外,不可能是其他有理数....超越数不能成为任何整数多项式的根..?
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN = 4e5 + 7;
const int P = 998244353;
const int G1 = 3;
const int G2 = (P + 1) / 3;
int N, len, L, R[MAXN];
ll F[MAXN], A[MAXN], B[MAXN], C[MAXN], D[MAXN];
ll fac[MAXN], ifac[MAXN], inv[MAXN];
inline ll ksm(ll x, ll y) {
ll ans = 1;
while(y) {
if(y & 1)ans = ans * x % P;
x = x * x % P;
y >>= 1;
}
return ans;
}
inline void init() {
fac[0] = 1;
for(int i = 1; i <= N; ++i)fac[i] = fac[i - 1] * i % P;
ifac[N] = ksm(fac[N], P - 2);
for(int i = N - 1; i >= 1; --i)ifac[i] = ifac[i + 1] * (i + 1) % P;
ifac[0] = 1;
for(int i = 1; i <= N; ++i)inv[i] = fac[i - 1] * ifac[i] % P;
return ;
}
inline void qiudao(ll *F, ll *B, int N) {
for(int i = 1; i < N; ++i) {
B[i - 1] = F[i] * i % P;
}
B[N - 1] = 0;
return ;
}
inline void jifen(ll *F, ll *B, int N) {
for(int i = N - 1; i >= 0; --i) {
B[i + 1] = F[i] * inv[i + 1] % P;
}
B[0] = 0;
return ;
}
inline void init(int x) {
len = 1;
L = 0;
while(len <= 2 * x)len <<= 1, ++L;
for(int i = 1; i < len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
}
inline void NTT(ll *F, int len, int typ) {
for(int i = 1; i < len; ++i)if(i < R[i])swap(F[i], F[R[i]]);
for(int mid = 1; mid < len; mid <<= 1) {
ll wn = ksm(typ == 1 ? G1 : G2, (P - 1) / (mid << 1));
for(int j = 0; j < len; j += (mid << 1)) {
ll w = 1;
for(int k = 0; k < mid; ++k, w = w * wn % P) {
int x = F[j + k], y = F[j + k + mid] * w % P;
F[j + k] = (x + y) % P;
F[j + k + mid] = (x + P - y) % P;
}
}
}
if(typ == -1) {
ll inv = ksm(len, P - 2);
for(int i = 0; i < len; ++i)F[i] = F[i] * inv % P;
}
return ;
}
inline void getinv(ll *F, ll *G, int n) {
if(n == 1) {
G[0] = ksm(F[0], P - 2);
return ;
}
getinv(F, G, (n + 1) / 2);
init(n);
for(int i = 0; i < n; ++i)D[i] = F[i], C[i] = G[i];
for(int i = n; i < len; ++i)D[i] = 0, C[i] = 0;
NTT(D, len, 1);
NTT(C, len, 1);
for(int i = 0; i < len; ++i)C[i] = (2 * C[i] % P - C[i] * C[i] % P * D[i] % P + P) % P;
NTT(C, len, -1);
for(int i = 0; i < n; ++i)G[i] = C[i];
return;
}
inline void getln(ll *F, int n) {
qiudao(F, A, n);
getinv(F, B, n);
init(n);
NTT(A, len, 1);
NTT(B, len, 1);
for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
NTT(A, len, -1);
jifen(A, F, n);
return ;
}
int main() {
scanf("%d", &N);
for(int i = 0; i < N; ++i)scanf("%lld", &F[i]);
init();
getln(F, N);
for(int i = 0; i < N; ++i)printf("%lld ", F[i]);
puts("");
return 0;
}
代码分析:没什么好说的,有一个很坑我的点,就是变量重名!
比如我们再ln中递归这个inv用A数组,然后inv的辅助数组又用了A....
其实解决方法就是多开两个数组,但是调试的时候可能很痛苦!
注意求导和积分顺序?
【模板】多项式指数函数(多项式 exp)
证明:
先同时取对数吧!
这个式子直接上牛顿迭代其实就做完了,具体的
牛顿迭代的式子是:
几何意义相当于我们先求出号点的切线方程,有
然后就是这个式子的根,带入y=0解得上式qwq
这个方法的精度增长的很快,基本上是10的次方增长的,对于迭代多项式,也是每次多两倍的精度
注意!我们这个F函数参数是一个多项式B,等价于A是一个常数,没有用!
因为的取值不随发生变化啊
下面我们有
于是我们对着最后一个式子做即可,时间复杂度就是的,当然,和lct一样的
注意,题目保证了如果没有这个条件是不能exp的,e^0才是1,其他都是超越数,常数项要是有理数
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int P = 998244353;
const int G1 = 3;
const int G2 = (P + 1) / 3;
const int MAXN = 4e5 + 7;
int N, len, L, R[MAXN];
ll F[MAXN], G[MAXN], A[MAXN], B[MAXN], C[MAXN], D[MAXN], E[MAXN], H[MAXN];
inline ll ksm(ll x, ll y) {
ll ans = 1;
while(y) {
if(y & 1)ans = ans * x % P;
x = x * x % P;
y >>= 1;
}
return ans;
}
ll ifac[MAXN], fac[MAXN], inv[MAXN];
inline void init2() {
fac[0] = 1;
for(int i = 1; i <= N; ++i)fac[i] = fac[i - 1] * i % P;
ifac[N] = ksm(fac[N], P - 2);
for(int i = N - 1; i >= 1; --i)ifac[i] = ifac[i + 1] * (i + 1) % P;
ifac[0] = 1;
for(int i = 1; i <= N; ++i)inv[i] = ifac[i] * fac[i - 1] % P;
return ;
}
inline void init(int x) {
L = 0;
len = 1;
while(len < 2 * x)len <<= 1, ++L;
for(int i = 0; i < len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
return ;
}
inline void NTT(ll *F, int len, int typ) {
for(int i = 0; i < len; ++i)if(i < R[i])swap(F[i], F[R[i]]);
for(int mid = 1; mid < len; mid <<= 1) {
ll wn = ksm(typ == 1 ? G1 : G2, (P - 1) / (mid << 1));
for(int j = 0; j < len; j += (mid << 1)) {
ll w = 1;
for(int k = 0; k < mid; ++k, w = w * wn % P) {
ll x = F[j + k], y = F[j + k + mid] * w % P;
F[j + k] = (x + y) % P;
F[j + k + mid] = (x + P - y) % P;
}
}
}
if(typ == -1) {
ll inv = ksm(len, P - 2);
for(int i = 0; i < len; ++i)F[i] = F[i] * inv % P;
}
return ;
}
//(G(x)-H(x))^2=0
//G(x)^2-2GH+H^2=0
//G(x)=2H-H^2F
inline void getinv(ll *F, ll *G, int n) {
if(n == 1) {
G[0] = ksm(F[0], P - 2);
return ;
}
getinv(F, G, (n + 1) / 2);
init(n);
for(int i = 0; i < n; ++i)A[i] = G[i], B[i] = F[i];
for(int i = n; i < len; ++i)A[i] = B[i] = 0;
NTT(A, len, 1);
NTT(B, len, 1);
for(int i = 0; i < len; ++i)A[i] = (2 * A[i] % P - A[i] * B[i] % P * A[i] % P + P) % P;
NTT(A, len, -1);
for(int i = 0; i < n; ++i)G[i] = A[i];
return ;
}
inline void jifen(ll *F, ll *B, int n) {
for(int i = n - 1; i >= 0; --i)B[i + 1] = F[i] * inv[i + 1] % P;
B[0] = 0;
return ;
}
inline void qiudao(ll *F, ll *B, int n) {
for(int i = 1; i < n; ++i) B[i - 1] = F[i] * i % P;
B[n - 1] = 0;
return ;
}
//(ln F(x))' = F'(x)/F(x)
inline void getln(ll *F, ll *G, int n) {
for(int i = 0; i <= n * 2; ++i)D[i] = C[i] = 0;//这句话是关键,如果我们那里init里写的是<=,这个就是错的
qiudao(F, C, n);
getinv(F, D, n);
init(n);
NTT(C, len, 1);
NTT(D, len, 1);
for(int i = 0; i < len; ++i)C[i] = C[i] * D[i] % P;
NTT(C, len, -1);
jifen(C, G, n);
return ;
}
//G_1(x)=G_0(x)*(1-ln (G_0(x))+F(x))) (mod x^n)
inline void getexp(ll *F, ll *G, int n) {
if(n == 1) {
G[0] = 1;
return ;
}
getexp(F, G, (n + 1) / 2);
getln(G, E, n);
init(n);
E[0] = (1 - E[0] + F[0] + P) % P;
for(int i = 1; i < n; ++i)E[i] = (-E[i] + F[i] + P) % P;
for(int i = 0; i < n; ++i)H[i] = G[i];
for(int i = n; i < len; ++i)H[i] = E[i] = 0;
NTT(H, len, 1);
NTT(E, len, 1);
for(int i = 0; i < len; ++i)E[i] = E[i] * H[i] % P;
NTT(E, len, -1);
for(int i = 0; i < n; ++i)G[i] = E[i];
for(int i = n; i < len; ++i)E[i] = G[i] = 0;
return ;
}
int main() {
scanf("%d", &N);
for(int i = 0; i < N; ++i)scanf("%lld", &F[i]);
init2();
getexp(F, G, N); //坏了!qwq
for(int i = 0; i < N; ++i) {
printf("%lld ", G[i]);
}
return 0;
}
代码分析:md一开始wqy他式子是个错的所以说我直接错到天上去了....
但是吴清月本人这个代码并没有问题,也就是说他的式子可能是等价的?
问题主要在于对于f'的求导,我们是对于f'而不是求导!
后来才检查出来一个巨大的锅:+1是常数项,不能每一项都+1
同时对于getln函数有个巨大的要求就是每次调用时要清空辅助数组,而且要清空到2倍,同时我们里面init数组的写法就要是那样子的就是<号而不是<=号,否则会getinv的时候循环卷积出锅再就是求导的时候n-1项清空而不是第n项/xk
以后都写<号吧,也不会有错,还能节省一点常数