与数论有关的模板
有一些数论相关的模板，经常是需要用到的。
下降幂多项式乘法
设
$$f(x)=\sum_{i=0}^na_ix^{\underline{i}}$$
$$G(x)=\sum_{i=0}^na_i\frac{x^i}{i!}$$
$$F(x)=\sum_{i=0}^nf(i)\frac{x^i}{i!}$$
其中$x^{\underline{n}}$表示$\prod_{i=0}^{n-1}(x-i)$，类似的还有$x^{\overline{n}}$表示上升幂$\prod_{i=0}^{n-1}(x+i)$。
即$A(x)$为原多项式，$G(x)$为系数的EGF，$F(x)$为点值的EGF。
我们发现$x^{\underline{n}}$对$F(x)$的贡献为$\sum_{i=n}^{\infty}\frac{i^{\underline{n}}}{i!}x^i=\sum_{i=n}^{\infty}\frac{1}{(i-n)!}x^i=x^ne^x$
所以$F=G*e^x$，即$G=F*e^{-x}$
点值就可以直接相乘了。
1 #include<bits/stdc++.h> 2 #define Rint register int 3 using namespace std; 4 typedef long long LL; 5 const int mod = 998244353, G = 3, Gi = 332748118, N = 800003; 6 inline int kasumi(int a, int b){ 7 int res = 1; 8 while(b){ 9 if(b & 1) res = (LL) res * a % mod; 10 a = (LL) a * a % mod; 11 b >>= 1; 12 } 13 return res; 14 } 15 int fac[N], inv[N]; 16 inline void init(int n){ 17 fac[0] = 1; 18 for(Rint i = 1;i <= n;i ++) fac[i] = (LL) i * fac[i - 1] % mod; 19 inv[n] = kasumi(fac[n], mod - 2); 20 for(Rint i = n;i;i --) inv[i - 1] = (LL) i * inv[i] % mod; 21 } 22 int rev[N]; 23 inline int calrev(int n){ 24 int limit = 1, L = -1; 25 while(limit <= n){limit <<= 1; L ++;} 26 for(Rint i = 0;i < limit;i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L); 27 return limit; 28 } 29 inline void NTT(int *A, int limit, int type){ 30 for(Rint i = 0;i < limit;i ++) if(i < rev[i]) swap(A[i], A[rev[i]]); 31 for(Rint mid = 1;mid < limit;mid <<= 1){ 32 int Wn = kasumi(type == 1 ? G : Gi, (mod - 1) / (mid << 1)); 33 for(Rint j = 0;j < limit;j += (mid << 1)){ 34 int w = 1; 35 for(Rint k = 0;k < mid;k ++, w = (LL) w * Wn % mod){ 36 int x = A[j + k], y = (LL) w * A[j + k + mid] % mod; 37 A[j + k] = (x + y) % mod; 38 A[j + k + mid] = (x - y + mod) % mod; 39 } 40 } 41 } 42 if(type == -1){ 43 int inv = kasumi(limit, mod - 2); 44 for(Rint i = 0;i < limit;i ++) A[i] = (LL) A[i] * inv % mod; 45 } 46 } 47 int n, m, A[N], B[N], C[N]; 48 int main(){ 49 scanf("%d%d", &n, &m); 50 for(Rint i = 0;i <= n;i ++) scanf("%d", A + i); 51 for(Rint i = 0;i <= m;i ++) scanf("%d", B + i); 52 init(n + m); 53 for(Rint i = 0;i <= n + m;i ++) C[i] = inv[i]; 54 int limit = calrev((n + m) << 1); 55 NTT(A, limit, 1); NTT(B, limit, 1); NTT(C, limit, 1); 56 for(Rint i = 0;i < limit;i ++){ 57 A[i] = (LL) A[i] * C[i] % mod; 58 B[i] = (LL) B[i] * C[i] % mod; 59 } 60 NTT(A, limit, -1); NTT(B, limit, -1); 61 for(Rint i = 0;i < limit;i ++) A[i] = (LL) A[i] * B[i] % mod * fac[i] % mod; 62 memset(C, 0, sizeof C); 63 for(Rint i = 0;i <= n + m;i ++) C[i] = (i & 1) ? (mod - inv[i]) : inv[i]; 64 NTT(A, limit, 1); NTT(C, limit, 1); 65 for(Rint i = 0;i < limit;i ++) A[i] = (LL) A[i] * C[i] % mod; 66 NTT(A, limit, -1); 67 for(Rint i = 0;i <= n + m;i ++) printf("%d ", A[i]); 68 }
点值转下降幂多项式
通常幂，下降幂，上升幂互相转化
通常幂转化为其他两个用第二类斯特林数.
转化为通常幂用第一类斯特林数。
小转大的时候带上$(-1)^{n-i}$的系数。
第二类斯特林数-行
$$\begin{Bmatrix}n \\ m\end{Bmatrix}=\sum_{i=0}^m\frac{(-1)^i}{i!}\times \frac{(m-i)^n}{(m-i)!}$$
1 #include<bits/stdc++.h> 2 #define Rint register int 3 using namespace std; 4 typedef long long LL; 5 const int N = 800003, mod = 167772161, g = 3, gi = 55924054; 6 int n, inv[N], A[N], B[N]; 7 inline int kasumi(int a, int b){ 8 int res = 1; 9 while(b){ 10 if(b & 1) res = (LL) res * a % mod; 11 a = (LL) a * a % mod; 12 b >>= 1; 13 } 14 return res; 15 } 16 int rev[N]; 17 inline void NTT(int *A, int limit, int type){ 18 for(Rint i = 0;i < limit;i ++) 19 if(i < rev[i]) swap(A[i], A[rev[i]]); 20 for(Rint mid = 1;mid < limit;mid <<= 1){ 21 int Wn = kasumi(type == 1 ? g : gi, (mod - 1) / (mid << 1)); 22 for(Rint j = 0;j < limit;j += mid << 1){ 23 int w = 1; 24 for(Rint k = 0;k < mid;k ++, w = (LL) w * Wn % mod){ 25 int x = A[j + k], y = (LL) w * A[j + k + mid] % mod; 26 A[j + k] = (x + y) % mod; 27 A[j + k + mid] = (x - y + mod) % mod; 28 } 29 } 30 } 31 if(type == -1){ 32 int inv = kasumi(limit, mod - 2); 33 for(Rint i = 0;i < limit;i ++) 34 A[i] = (LL) A[i] * inv % mod; 35 } 36 } 37 int main(){ 38 scanf("%d", &n); 39 inv[1] = 1; 40 for(Rint i = 2;i <= n;i ++) inv[i] = (LL) (mod - mod / i) * inv[mod % i] % mod; 41 inv[0] = 1; 42 for(Rint i = 1;i <= n;i ++) inv[i] = (LL) inv[i] * inv[i - 1] % mod; 43 for(Rint i = 0;i <= n;i ++){ 44 A[i] = (i & 1) ? -inv[i] : inv[i]; 45 B[i] = (LL) kasumi(i, n) * inv[i] % mod; 46 } 47 int limit = 1, L = -1; 48 while(limit <= (n + 1 << 1)){limit <<= 1; L ++;} 49 for(Rint i = 0;i < limit;i ++) 50 rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L); 51 NTT(A, limit, 1); NTT(B, limit, 1); 52 for(Rint i = 0;i < limit;i ++) 53 A[i] = (LL) A[i] * B[i] % mod; 54 NTT(A, limit, -1); 55 for(Rint i = 0;i <= n;i ++) 56 printf("%d ", A[i]); 57 }
第二类斯特林数-列
$$S_m(x)=\sum_{i=0}^{+\infty}x^i\begin{Bmatrix}i \\ m\end{Bmatrix}$$
$$\begin{Bmatrix}n \\ m\end{Bmatrix}=\begin{Bmatrix}n-1 \\ m-1\end{Bmatrix}+\begin{Bmatrix}n-1 \\ m\end{Bmatrix}\times m$$
$$S_m(x)=mxS_m(x)+xS_{m-1}(x)$$
$$S_m(x)=S_{m-1}(x)\times \frac{x}{1-mx}$$
$$S_m(x)=\frac{x^m}{\prod_{i=1}^m(1-ix)}$$
分治计算连乘积，$O(n\log^2 n)$。
$$\sum_{i=k}^{+\infty}\begin{Bmatrix}i \\ k\end{Bmatrix}\frac{x^i}{i!}=\frac{(e^x-1)^k}{k!}$$
这个式子的组合意义就是，左边是把$i$个无区别的元素分成$k$个无序的集合，右边表示每个集合的EGF是$e^x-1$，然后去标号是除以$k!$。
虽然这个做法是$O(n\log n)$，但是实际上更慢。（exp大常数）
1 #include<bits/stdc++.h> 2 #define Rint register int 3 using namespace std; 4 typedef long long LL; 5 const int N = 1 << 19, mod = 167772161, g = 3, gi = 55924054; 6 int n, m; 7 inline int kasumi(int a, int b){ 8 int res = 1; 9 while(b){ 10 if(b & 1) res = (LL) res * a % mod; 11 a = (LL) a * a % mod; 12 b >>= 1; 13 } 14 return res; 15 } 16 int rev[N]; 17 inline void NTT(int *A, int limit, int type){ 18 for(Rint i = 0;i < limit;i ++) 19 if(i < rev[i]) swap(A[i], A[rev[i]]); 20 for(Rint mid = 1;mid < limit;mid <<= 1){ 21 int Wn = kasumi(type == 1 ? g : gi, (mod - 1) / (mid << 1)); 22 for(Rint j = 0;j < limit;j += (mid << 1)){ 23 int w = 1; 24 for(Rint k = 0;k < mid;k ++, w = (LL) w * Wn % mod){ 25 int x = A[j + k], y = (LL) A[j + k + mid] * w % mod; 26 A[j + k] = (x + y) % mod; 27 A[j + k + mid] = (x - y + mod) % mod; 28 } 29 } 30 } 31 if(type == -1){ 32 int inv = kasumi(limit, mod - 2); 33 for(Rint i = 0;i < limit;i ++) A[i] = (LL) A[i] * inv % mod; 34 } 35 } 36 int ans[N]; 37 inline void poly_inv(int *A, int deg){ 38 static int tmp[N]; 39 if(deg == 1){ 40 ans[0] = kasumi(A[0], mod - 2); 41 return; 42 } 43 poly_inv(A, deg + 1 >> 1); 44 int limit = 1, L = -1; 45 while(limit <= (deg << 1)){limit <<= 1; L ++;} 46 for(Rint i = 0;i < limit;i ++) 47 rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L); 48 for(Rint i = 0;i < deg;i ++) tmp[i] = A[i]; 49 for(Rint i = deg;i < limit;i ++) tmp[i] = 0; 50 NTT(tmp, limit, 1); NTT(ans, limit, 1); 51 for(Rint i = 0;i < limit;i ++) 52 ans[i] = (2 - (LL) tmp[i] * ans[i] % mod + mod) % mod * ans[i] % mod; 53 NTT(ans, limit, -1); 54 for(Rint i = deg;i < limit;i ++) ans[i] = 0; 55 } 56 int A[20][N]; 57 inline void calc(int dep, int L, int R){ 58 if(L == R){ 59 A[dep][0] = 1; A[dep][1] = mod - L; 60 return; 61 } 62 int mid = L + R >> 1; 63 calc(dep + 1, L, mid); 64 for(Rint i = 0;i <= mid - L + 1;i ++) A[dep][i] = A[dep + 1][i]; 65 calc(dep + 1, mid + 1, R); 66 int limit = 1, l = -1; 67 while(limit <= (R - L + 1)){limit <<= 1; l ++;} 68 for(Rint i = 0;i < limit;i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l); 69 for(Rint i = mid - L + 2;i < limit;i ++) A[dep][i] = 0; 70 for(Rint i = R - mid + 1;i < limit;i ++) A[dep + 1][i] = 0; 71 NTT(A[dep], limit, 1); NTT(A[dep + 1], limit, 1); 72 for(Rint i = 0;i < limit;i ++) A[dep][i] = (LL) A[dep][i] * A[dep + 1][i] % mod; 73 NTT(A[dep], limit, -1); 74 } 75 int main(){ 76 // freopen("my.out", "w", stdout); 77 scanf("%d%d", &n, &m); 78 if(n < m){ 79 for(Rint i = 0;i <= n;i ++) printf("0 "); 80 return 0; 81 } 82 calc(0, 1, m); 83 // for(Rint i = 0;i <= m;i ++) printf("%d ", A[0][i]); 84 // putchar('\n'); 85 poly_inv(A[0], n - m + 1); 86 for(Rint i = 0;i < m;i ++) printf("0 "); 87 for(Rint i = 0;i <= n - m;i ++) printf("%d ", ans[i]); 88 }
第一类斯特林数-列
$$(x+1)^k=\sum_{i=0}^kx^i\frac{k^{\underline{i}}}{i!}$$
$$=\sum_{i=0}^k\frac{x^i}{i!}\sum_{j=0}^i(-1)^{i-j}\begin{bmatrix}i \\ j\end{bmatrix}k^j$$
$$=\sum_{i=0}^nk^i\sum_{j=i}^k(-1)^{j-i}\begin{bmatrix}j \\ i\end{bmatrix}\frac{x^j}{j!}$$
用另一种方法展开可以得到：
$$(x+1)^k=\exp(k\ln(x+1))=\sum_{i=0}^{+\infty}k^i\frac{\ln(x+1)^i}{i!}$$
所以
$$\sum_{i=k}^\infty(-1)^{i-k}\begin{bmatrix}i \\ k\end{bmatrix}\frac{x^i}{i!}=\frac{\ln(x+1)^k}{k!}$$
这里先贴一个90分代码，之后再卡常。。。
1 #include<bits/stdc++.h> 2 #define Rint register int 3 using namespace std; 4 typedef long long LL; 5 const int mod = 167772161, G = 3, Gi = 55924054, N = 1 << 19; 6 inline int kasumi(int a, int b){ 7 int res = 1; 8 while(b){ 9 if(b & 1) res = (LL) res * a % mod; 10 a = (LL) a * a % mod; 11 b >>= 1; 12 } 13 return res; 14 } 15 int fac[N], inv[N], invfac[N]; 16 inline void init(int n){ 17 fac[0] = 1; 18 for(Rint i = 1;i <= n;i ++) fac[i] = (LL) i * fac[i - 1] % mod; 19 invfac[n] = kasumi(fac[n], mod - 2); 20 for(Rint i = n;i;i --){ 21 invfac[i - 1] = (LL) i * invfac[i] % mod; 22 inv[i] = (LL) invfac[i] * fac[i - 1] % mod; 23 } 24 } 25 int rev[N]; 26 inline int calrev(int n){ 27 int limit = 1, L = -1; 28 while(limit <= n){limit <<= 1; L ++;} 29 for(Rint i = 0;i < limit;i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L); 30 return limit; 31 } 32 inline void NTT(int *A, int limit, int type){ 33 for(Rint i = 0;i < limit;i ++) if(i < rev[i]) swap(A[i], A[rev[i]]); 34 for(Rint mid = 1;mid < limit;mid <<= 1){ 35 int Wn = kasumi(type == 1 ? G : Gi, (mod - 1) / (mid << 1)); 36 for(Rint j = 0;j < limit;j += (mid << 1)){ 37 int w = 1; 38 for(Rint k = 0;k < mid;k ++, w = (LL) w * Wn % mod){ 39 int x = A[j + k], y = (LL) w * A[j + k + mid] % mod; 40 A[j + k] = (x + y) % mod; 41 A[j + k + mid] = (x - y + mod) % mod; 42 } 43 } 44 } 45 if(type == -1){ 46 int inv = kasumi(limit, mod - 2); 47 for(Rint i = 0;i < limit;i ++) A[i] = (LL) A[i] * inv % mod; 48 } 49 } 50 int ans[N]; 51 inline void poly_inv(int *A, int deg){ 52 static int tmp[N]; 53 if(deg == 1){ 54 ans[0] = kasumi(A[0], mod - 2); 55 return; 56 } 57 poly_inv(A, deg + 1 >> 1); 58 int limit = calrev(deg << 1); 59 for(Rint i = 0;i < deg;i ++) tmp[i] = A[i]; 60 for(Rint i = deg;i < limit;i ++) tmp[i] = ans[i] = 0; 61 NTT(tmp, limit, 1); NTT(ans, limit, 1); 62 for(Rint i = 0;i < limit;i ++) ans[i] = (2 - (LL) ans[i] * tmp[i] % mod + mod) * ans[i] % mod; 63 NTT(ans, limit, -1); 64 for(Rint i = deg;i < limit;i ++) ans[i] = 0; 65 } 66 int Ln[N]; 67 inline void poly_Ln(int *A, int deg){ 68 static int tmp[N]; 69 if(deg == 1){ 70 Ln[0] = 0; 71 return; 72 } 73 poly_inv(A, deg); 74 int limit = calrev(deg << 1); 75 for(Rint i = 1;i < deg;i ++) tmp[i - 1] = (LL) i * A[i] % mod; 76 for(Rint i = deg - 1;i < limit;i ++) tmp[i] = 0; 77 NTT(ans, limit, 1); NTT(tmp, limit, 1); 78 for(Rint i = 0;i < limit;i ++) Ln[i] = (LL) ans[i] * tmp[i] % mod; 79 NTT(Ln, limit, -1); 80 for(Rint i = deg + 1;i < limit;i ++) Ln[i] = 0; 81 for(Rint i = deg;i;i --) Ln[i] = (LL) Ln[i - 1] * inv[i] % mod; 82 Ln[0] = 0; 83 for(Rint i = 0;i < limit;i ++) ans[i] = tmp[i] = 0; 84 } 85 int Exp[N]; 86 inline void poly_Exp(int *A, int deg){ 87 if(deg == 1){ 88 Exp[0] = 1; 89 return; 90 } 91 poly_Exp(A, deg + 1 >> 1); 92 poly_Ln(Exp, deg); 93 for(Rint i = 0;i < deg;i ++) Ln[i] = (A[i] + (i == 0) - Ln[i] + mod) % mod; 94 int limit = calrev(deg << 1); 95 NTT(Exp, limit, 1); NTT(Ln, limit, 1); 96 for(Rint i = 0;i < limit;i ++) Exp[i] = (LL) Exp[i] * Ln[i] % mod; 97 NTT(Exp, limit, -1); 98 for(Rint i = deg;i < limit;i ++) Exp[i] = ans[i] = Ln[i] = 0; 99 for(Rint i = 0;i < deg;i ++) ans[i] = Ln[i] = 0; 100 } 101 int n, k, A[N]; 102 int main(){ 103 scanf("%d%d", &n, &k); 104 if(n < k){ 105 for(Rint i = 0;i <= n;i ++) printf("0 "); 106 return 0; 107 } 108 init(n); 109 A[0] = A[1] = 1; 110 poly_Ln(A, n + 1); 111 // for(Rint i = 0;i <= n;i ++) printf("%d ", Ln[i]); 112 // putchar('\n'); 113 int lst = -1, val = 0; 114 for(Rint i = 0;i <= n;i ++) if(Ln[i]) {val = Ln[lst = i]; break;} 115 for(Rint i = 0;i <= n - lst;i ++) A[i] = (LL) Ln[i + lst] * inv[val] % mod; 116 poly_Ln(A, n + 1); 117 for(Rint i = 0;i <= n;i ++) A[i] = (LL) k * Ln[i] % mod; 118 poly_Exp(A, n + 1); 119 val = kasumi(val, k); 120 for(Rint i = 0;i < lst * k;i ++) A[i] = 0; 121 for(Rint i = 0;i <= n - lst * k;i ++) A[i + lst * k] = (LL) val * Exp[i] % mod * invfac[k] % mod; 122 for(Rint i = 0;i <= n;i ++) 123 printf("%d ", (LL) A[i] * fac[i] % mod * ((i - k) & 1 ? (mod - 1) : 1) % mod); 124 }
第一类斯特林数-行
$$\sum_{i=0}^{+\infty}\begin{bmatrix}n \\ i\end{bmatrix}x^i=\prod_{i=0}^{n-1}(x+i)$$
1 #include<bits/stdc++.h> 2 #define Rint register int 3 using namespace std; 4 typedef long long LL; 5 const int N = 1 << 19, mod = 167772161, g = 3, gi = 55924054; 6 inline int kasumi(int a, int b){ 7 int res = 1; 8 while(b){ 9 if(b & 1) res = (LL) res * a % mod; 10 a = (LL) a * a % mod; 11 b >>= 1; 12 } 13 return res; 14 } 15 int fac[N], inv[N], invfac[N]; 16 inline void init(int n){ 17 fac[0] = 1; 18 for(Rint i = 1;i <= n;i ++) fac[i] = (LL) i * fac[i - 1] % mod; 19 invfac[n] = kasumi(fac[n], mod - 2); 20 for(Rint i = n;i;i --){ 21 invfac[i - 1] = (LL) i * invfac[i] % mod; 22 inv[i] = (LL) fac[i - 1] * invfac[i] % mod; 23 } 24 } 25 int rev[N]; 26 inline int calrev(int n){ 27 int limit = 1, L = -1; 28 while(limit <= n){limit <<= 1; L ++;} 29 for(Rint i = 0;i < limit;i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L); 30 return limit; 31 } 32 inline void NTT(int *A, int limit, int type){ 33 for(Rint i = 0;i < limit;i ++) if(i < rev[i]) swap(A[i], A[rev[i]]); 34 for(Rint mid = 1;mid < limit;mid <<= 1){ 35 int Wn = kasumi(type == 1 ? g : gi, (mod - 1) / (mid << 1)); 36 for(Rint j = 0;j < limit;j += (mid << 1)){ 37 int w = 1; 38 for(Rint k = 0;k < mid;k ++, w = (LL) w * Wn % mod){ 39 int x = A[j + k], y = (LL) w * A[j + k + mid] % mod; 40 A[j + k] = (x + y) % mod; 41 A[j + k + mid] = (x - y + mod) % mod; 42 } 43 } 44 } 45 if(type == -1){ 46 int inv = kasumi(limit, mod - 2); 47 for(Rint i = 0;i < limit;i ++) A[i] = (LL) A[i] * inv % mod; 48 } 49 } 50 int F[N], G[N], tmp[N]; 51 inline void work(int n){ 52 if(n == 1){ 53 F[0] = 0; F[1] = 1; 54 return; 55 } 56 int mid = n >> 1; work(mid); 57 for(Rint i = 0, now = 1;i <= mid;i ++, now = (LL) now * mid % mod){ 58 tmp[i] = (LL) F[i] * fac[i] % mod; 59 G[mid - i] = (LL) now * invfac[i] % mod; 60 } 61 int limit = calrev(n); 62 NTT(tmp, limit, 1); NTT(G, limit, 1); 63 for(Rint i = 0;i < limit;i ++) tmp[i] = (LL) G[i] * tmp[i] % mod; 64 NTT(tmp, limit, -1); 65 for(Rint i = 0;i <= mid;i ++) 66 G[i] = (LL) tmp[mid + i] * invfac[i] % mod; 67 for(Rint i = mid + 1;i < limit;i ++) F[i] = G[i] = 0; 68 NTT(F, limit, 1); NTT(G, limit, 1); 69 for(Rint i = 0;i < limit;i ++) F[i] = (LL) F[i] * G[i] % mod; 70 NTT(F, limit, -1); 71 for(Rint i = (mid << 1) + 1;i < limit;i ++) F[i] = 0; 72 for(Rint i = 0;i < limit;i ++) G[i] = tmp[i] = 0; 73 if(n & 1){ 74 for(Rint i = 0;i <= n;i ++) tmp[i] = (LL) (n - 1) * F[i] % mod; 75 for(Rint i = 0;i < n;i ++) tmp[i + 1] = (tmp[i + 1] + F[i]) % mod; 76 for(Rint i = 0;i <= n;i ++) F[i] = tmp[i]; 77 } 78 } 79 int n; 80 int main(){ 81 scanf("%d", &n); 82 init(n); work(n); 83 for(Rint i = 0;i <= n;i ++) printf("%d ", F[i]); 84 }
