devbox
爱喝兽奶帝天荒
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

【学习笔记】多项式全家桶（包含全套证明）_多项式全家桶csdn

作者：爱喝兽奶帝天荒 | 2024-08-08 12:47:38

赞

踩

多项式全家桶csdn

整理的算法模板合集： ACM模板

点我看算法全家桶系列！！！

实际上是一个全新的精炼模板整合计划

多项式全家桶！！！

目录

在这里插入图片描述
tips：

注意一点，在我们每次求 limit 的时候，把范围都乘 2 ，反正乘了开大一点不会错，最多会跑的慢一点，但是有时候不开就会 WA（比如在求多项式除法的时候）
limit 可以设成全局变量，注意每次使用 limit 的时候都要按照下面的格式初始化一下 limit 以及 L
应该是NTT
就是这里：

       for(limit = 1, L = 0; limit <= (n + m) * 2; limit <<= 1) L ++ ;
       for(int i = 0; i < limit; ++ i) 
           RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));

  • 1
  • 2
  • 3

0x00.多项式基本概述

既然是多项式全家桶，那么肯定要先来了解一下什么是多项式hhh，基础一定要扎实，不然到后面让你怀疑人生

0x01. 多项式的基础概念

首先是多项式的一些基础的概念：
学习笔记

0x02. 多项式的度

对于一个多项式 f ( x ) f(x) f(x)，称其最高次项的次数为该多项式的 度（Degree） ，记作 d e g   f deg\ f deg f。

好吧就是多项式的次数。（顺带一提，单项式的次数是所有字母次数之和hhh）

0x03. 多项式的乘法

我们上面也说了多项式 × \times × 多项式 = = = 新的多项式。

所以我们最关心的操作就是两个多项式的乘法 —— 卷积。

广义上的卷积为：

h ( x ) = ∫ − ∞ ∞ g ( τ ) ⋅ f ( x − τ ) d τ {h(x) = \int _{- \infty} ^{\infty}g(\tau) \cdot f(x - \tau)} \rm{d\tau} h(x)=∫−∞∞​g(τ)⋅f(x−τ)dτ

我们这里仅讨论多项式域，即两个多项式 / 数组 / 序列的卷积。

即给定两个多项式 f ( x ) f(x) f(x)， g ( x ) g(x) g(x) ：

f ( x ) = a 0 + a 1 x + ⋯ + a n x n f(x)= a_0 + a_1x+ \dots +a_nx^n f(x)=a0​+a1​x+⋯+an​xn

g ( x ) = b 0 + b 1 x + ⋯ + b m x m g(x)= b_0 + b_1x+ \dots +b_mx^m g(x)=b0​+b1​x+⋯+bm​xm

计算多项式 Q ( x ) = f ( x )   ⋅   g ( x ) : Q(x) = f(x)\ ·\ g(x): Q(x)=f(x) ⋅ g(x):

Q ( x ) = ∑ i = 0 n ∑ j = 0 m a i b j x i + j = c 0 + c 1 x + ⋯ + c n + m x n + m \boxed {Q(x) = \sum \limits_ {i = 0} ^ n \sum \limits_ {j = 0 } ^ m a_i b_j x ^ {i + j}} = c_0 + c_1 x + \dots + c_ {n + m} x ^ {n + m} Q(x)=i=0∑n​j=0∑m​ai​bj​xi+j​=c0​+c1​x+⋯+cn+m​xn+m

C ( x ) = A ( x ) ∗ B ( x ) = ∑ i = 0 n ( ∑ j = 0 i a j b i − j ) x j + i − j = ∑ i = 0 n ( ∑ j = 0 i a j b i − j ) x i C(x)=A(x)*B(x)=\sum_{i = 0}^{n}(\sum_{j=0}^{i}a_j b_{i-j})x^{j+i-j}=\sum_{i = 0}^{n}(\sum_{j=0}^{i}a_j b_{i-j})x^{i} C(x)=A(x)∗B(x)=i=0∑n​(j=0∑i​aj​bi−j​)xj+i−j=i=0∑n​(j=0∑i​aj​bi−j​)xi

也就是对于 C ( x ) C(x) C(x) 的第 i i i 项的系数：

[ i ] C ( x ) = ∑ j = 0 i a j b i − j [i]C(x)=\sum_{j=0}^{i}a_j b_{i-j} [i]C(x)=j=0∑i​aj​bi−j​

C ( x ) = A ( x ) ∗ B ( x ) = ∑ k = 0 n + m − 2 ( ∑ k = i + j a i b j ) x k C(x)=A(x)* B(x)=\sum_{k=0}^{n+m-2}(\sum_{k=i+j}a_ib_j)x^k C(x)=A(x)∗B(x)=k=0∑n+m−2​(k=i+j∑​ai​bj​)xk

       这里使用的就是上面介绍的系数表示，我们如果直接暴力相乘很明显是一个 O ( n m ) O(nm) O(nm) 的复杂度，因为我们每一项都需要相乘。但是我们可以通过快速傅里叶变换将常用的系数表示转为点值表示用 O ( n l o g n ) O(nlogn) O(nlogn) 的时间复杂度下计算得到答案，再逆变换转成系数表示输出，将会在下面章节介绍。

0x04. 多项式的逆元

对于 n n n 次多项式 f ( x ) f(x) f(x)，若存在 g ( x ) g(x) g(x)，满足：
f ( x ) g ( x ) ≡ 1 ( m o d x n )   & &   deg ⁡ ⁡ g ≤ deg ⁡ ⁡ f f(x) g(x) \equiv 1 \pmod{x^{n}}\ \&\&\ \operatorname{\deg}{g} \le \operatorname{\deg}{f} f(x)g(x)≡1(modxn) && degg≤degf
则称 g ( x ) g(x) g(x)是 n n n 次多项式 f ( x ) f(x) f(x) 的逆元。

0x05. 多项式多点求值与插值

在这里插入图片描述

0x06 卷积的定义

对于一个序列，将其中元素一一映射到一个多项式函数的系数上， 这个多项式函数便叫做该序列的生成函数。

形式化地讲，对于序列 f 0 , f 1 , ⋯   , f n − 1 f_0,f_1,\cdots,f_{n-1} f0​,f1​,⋯,fn−1​ ， f ( x ) = ∑ k = 0 n − 1 f k x k f(x)=\displaystyle\sum_{k=0}^{n-1}f_kx^k f(x)=k=0∑n−1​fk​xk 为其生成函数。

卷积即为生成函数的乘积在对应序列的变换上的的抽象，“卷”即为其作用效果，“积”即为其本质。

对于序列 f , g f,g f,g ，其卷积序列 f ⊗ g f\otimes g f⊗g 满足 ( f ⊗ g ) k = ∑ i = 0 k f i × g k − i = ∑ i , j i + j = k f i × g j (f\otimes g)_k=\displaystyle\sum\limits_{i=0}^kf_i\times g_{k-i}=\sum\limits_{i,j}^{i+j=k}f_i\times g_j (f⊗g)k​=i=0∑k​fi​×gk−i​=i,j∑i+j=k​fi​×gj​

对于多项式 f , g f,g f,g，其多项式的卷积为： f ⊗ g = ∑ k = 0 n ( ∑ i , j i + j = k a i × b j ) x k f\otimes g=\sum\limits_{k=0}^{n}(\sum\limits_{i,j}^{i+j=k}a_i\times b_j)x^k f⊗g=k=0∑n​(i,j∑i+j=k​ai​×bj​)xk

而我们所熟知的 FFT 计算的是循环卷积，也即 ( f ⊗ g ) k = ∑ i + j ≡ k ( m o d n ) f i × g j (f\otimes g)_k=\displaystyle\sum_{i+j\equiv k\pmod n}f_i\times g_j (f⊗g)k​=i+j≡k(modn)∑​fi​×gj​ ， n n n 为序列长度。

0x07 卷积的基本性质

这里的卷积均为序列的卷积，因此证明时我们使用数学归纳法，即证明第 k k k 项成立，则整个序列均成立。由于多项式实际上就是序列的生成函数，所以性质同样成立。

  • f ⊗ g = g ⊗ f f\otimes g=g\otimes f f⊗g=g⊗f（交换律）

    我们使用定义 ( f ⊗ g ) k = ∑ i + j = k f i × g j (f\otimes g)_k=\displaystyle\sum_{i+j=k}f_i\times g_j (f⊗g)k​=i+j=k∑​fi​×gj​ ，以及乘法交换律 a × b = b × a a\times b=b\times a a×b=b×a 即可证明，因为 i + j = k i+j=k i+j=k 的 i i i 和 j j j 是一一对应的。

  • ( f ⊗ g ) ⊗ h = f ⊗ ( g ⊗ h ) (f\otimes g)\otimes h=f\otimes(g \otimes h) (f⊗g)⊗h=f⊗(g⊗h) （结合律）
    证明：
    [ f ⊗ ( g ⊗ h ) ] n = ∑ i = 0 n f n − i × ( g ⊗ h ) i = ∑ i = 0 n f n − i ∑ j = 0 i g j h i − j = ∑ i = 0 n ∑ j = 0 i f n − i g j h i − j = ∑ i + j + k = n f i g j h k [f⊗(g⊗h)]n=n∑i=0fn−i×(g⊗h)i=n∑i=0fn−ii∑j=0gjhi−j=n∑i=0i∑j=0fn−igjhi−j=∑i+j+k=nfigjhk [f⊗(g⊗h)]n​​=i=0∑n​fn−i​×(g⊗h)i​=i=0∑n​fn−i​j=0∑i​gj​hi−j​=i=0∑n​j=0∑i​fn−i​gj​hi−j​=i+j+k=n∑​fi​gj​hk​​

交换顺序后反向推导即可得证。

  • ( f ⊕ g ) ⊗ h = ( f ⊗ h ) ⊕ ( g ⊗ h ) (f\oplus g)\otimes h=(f\otimes h)\oplus(g\otimes h) (f⊕g)⊗h=(f⊗h)⊕(g⊗h) （分配律）
    其中 ( f ⊕ g ) k = a f k + b g k (f\oplus g)_k=af_k+bg_k (f⊕g)k​=afk​+bgk​ ，即序列 f , g f,g f,g 的加法， a , b a,b a,b 为常数。

证明：
[ ( f ⊕ g ) ⊗ h ] k = ∑ i = 0 k ( f ⊕ g ) i h k − i = ∑ i = 0 k ( a f i + b g i ) h k − i = ∑ i = 0 k ( a f i h k − i + b g i h k − i ) = a ∑ i = 0 k f i h k − i + b ∑ i = 0 k g i h k − i = a ( f ⊗ h ) k + b ( g ⊗ h ) k = [ ( f ⊗ h ) ⊕ ( g ⊗ h ) ] k [(f⊕g)⊗h]k=k∑i=0(f⊕g)ihk−i=k∑i=0(afi+bgi)hk−i=k∑i=0(afihk−i+bgihk−i)=ak∑i=0fihk−i+bk∑i=0gihk−i=a(f⊗h)k+b(g⊗h)k=[(f⊗h)⊕(g⊗h)]k [(f⊕g)⊗h]k​​=i=0∑k​(f⊕g)i​hk−i​=i=0∑k​(afi​+bgi​)hk−i​=i=0∑k​(afi​hk−i​+bgi​hk−i​)=ai=0∑k​fi​hk−i​+bi=0∑k​gi​hk−i​=a(f⊗h)k​+b(g⊗h)k​=[(f⊗h)⊕(g⊗h)]k​​

0x10.多项式插值

0x11.拉格朗日插值

0x11.1 拉格朗日插值简述

定理： n + 1 n+1 n+1 个点可以唯一确定一个 n n n 次多项式 f ( x ) f(x) f(x)。

现在给出n个点 ( x i , y i ) (x_i,y_i) (xi​,yi​)，请确定多项式 f ( k )   m o d   988244535 f(k)\ mod\ 988244535 f(k) mod 988244535的值。
在这里插入图片描述

0x11.2 拉格朗日插值模板

我们按照上面的思路，带入n个点求一下 f ( k ) f(k) f(k)即可。时间复杂度 O ( n 2 ) O(n^2) O(n2)。

注意本题还要求逆元，为了防止求逆元的时间复杂度影响整体的时间复杂度，所以我们分别计算出分子和分母，再将分子乘进分母的逆元，累加进最后的答案，时间复杂度的瓶颈就不会在求逆元上，总体的时间复杂度为 O ( n 2 ) O(n^2) O(n2)。

代码实现：

本题需要累乘，会爆int，记得开long long
（数学题开long long是个好习惯 ）

const int mod = 998244353;
ll n, m, k;
struct Point
{
    ll x, y;
}A[N];
ll inv(ll x) {return qpow(x, mod - 2, mod);}//快速幂的代码我就删了省空间，应该都会

int main()
{
    scanf("%lld%lld", &n, &k);
    for(int i = 1; i <= n; ++ i) {
        scanf("%lld%lld", &A[i].x, &A[i].y);
    }
    ll ans = 0;

    for(int i = 1; i <= n; ++ i) {
        ll s1 = A[i].y % mod;
        ll s2 = 1ll;
        for(int j = 1; j <= n; ++ j) {
            if(i != j) {
                s1 = s1 * (k - A[j].x) % mod;
                s2 = s2 * (A[i].x - A[j].x) % mod;
            }
        }
        ans += s1 * inv(s2) % mod;
    }
    printf("%lld\n", (ans % mod + mod) % mod);
    return 0;
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30

我们使用拉格朗日插值公式，对于一个数 k k k ，我们很容易在 O ( n 2 ) O(n^2) O(n2) 时间求得 F ( k ) F(k) F(k) 的数值，如果 x i x_i xi​ 是连续的，我们甚至可以利用预处理在 O ( n ) O(n) O(n) 时间内得到 F ( k ) F(k) F(k) 的数值。
但是如果 x i x_i xi​ 不连续，又有多组查询，就需要得到这个多项式的系数以保证求一个函数值的时间为 O ( n ) O(n) O(n) 。

0x11.3 拉格朗日插值法求系数

首先观察式子，发现对于每个 i i i , 上面部分是 ( x − x 1 ) × ( x − x 2 ) … × ( x − x n ) ( x − x i ) \cfrac{(x-x_1) \times (x-x_2)… \times (x-x_n)}{(x-x_i)} (x−xi​)(x−x1​)×(x−x2​)…×(x−xn​)​, 下面那部分和 y i y_i yi​ 都是常数。

所以可以 O ( n 2 ) O(n^2) O(n2) 处理出 ( x − x 1 ) × ( x − x 2 ) … × ( x − x n ) (x-x_1) \times (x-x_2)… \times (x-x_n) (x−x1​)×(x−x2​)…×(x−xn​) 这个 n n n 次多项式，然后通过模拟长除法， O ( n ) O(n) O(n)时间内可以得到 ( x − x 1 ) × ( x − x 2 ) … × ( x − x n ) ( x − x i ) \cfrac{(x-x_1) \times (x-x_2)… \times (x-x_n)}{(x-x_i)} (x−xi​)(x−x1​)×(x−x2​)…×(x−xn​)​，然后常系数直接 O ( n ) O(n) O(n) 暴力算出来，就得到了 x i x_i xi​ 对应的多项式，最后把所有多项式加起来就得到了最终的系数。

const int maxn = 5007;
ll mod = 998244353;

ll qpow(ll a, ll b) {
    ll res = 1;
    while (b) {
        if (b & 1)
            res = res * a % mod;
        a = a * a % mod, b >>= 1;
    }
    return res;
}
ll a[maxn], b[maxn], c[maxn], temp[maxn];
ll x[maxn], y[maxn];
int n;
void mul(ll *f, int len, ll t) { //len为多项式的次数+1，函数让多项式f变成f*(x+t)
    for (int i = len; i > 0; --i)
        temp[i] = f[i], f[i] = f[i - 1];
    temp[0] = f[0], f[0] = 0;
    for (int i = 0; i <= len; ++i)
        f[i] = (f[i] + t * temp[i]) % mod;
}
void dev(ll *f, ll *r, ll t) { //f是被除多项式的系数，r保存f除以x+t的结果
    for (int i = 0; i <= n; ++i)
        temp[i] = f[i];
    for (int i = n; i > 0; --i) {
        r[i - 1] = temp[i];
        temp[i - 1] = (temp[i - 1] - t * temp[i]) % mod;
    }
    return;
}
void lglr() {
    memset(a, 0, sizeof a);
    b[1] = 1, b[0] = -x[1];
    for (int i = 2; i <= n; ++i) {
        mul(b, i, -x[i]);
    }//预处理(x-x1)*(x-x2)...*(x-xn)
    for (int i = 1; i <= n; ++i) {
        ll fz = 1;
        for (int j = 1; j <= n; ++j) {
            if (j == i)
                continue;
            fz = fz * (x[i] - x[j]) % mod;
        }
        fz = qpow(fz, mod - 2);
        fz = fz * y[i] % mod; //得到多项式系数
        dev(b, c, -x[i]);//得到多项式，保存在b数组
        for (int j = 0; j < n; ++j)
            a[j] = (a[j] + fz * c[j]) % mod;
    }
}
int main() {
    ll k;
    cin >> n >> k;
    for (int i = 1; i <= n; ++i)
        scanf("%lld%lld", &x[i], &y[i]);
    lglr();
    ll ans = 0;
    ll res = 1;
    for (int i = 0; i < n; ++i) {
        ans = (ans + res * a[i]) % mod;
        res = res * k % mod;
    }
    ans = (ans + mod) % mod;
    cout << ans << endl;
}


  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67

拉格朗日插值法 在ACM竞赛中的一些应用

0x12.牛顿插值

待更hhh

关于多项式的一点研究及其在ACM竞赛中的应用

0x20.多项式乘法

0x21.快速傅里叶变换（FFT）

如果还没有学FFT的话可以点开下面链接学习一下，上万字长文带你速通FFT

0x21.1 FFT 简述

我们已知一个 ( n − 1 ) (n - 1) (n−1) 次多项式 A ( x )   =   ∑ i = 0 n − 1 a i x i A(x)~=~\sum_{i = 0}^{n - 1} a_i x^i A(x) = ∑i=0n−1​ai​xi 进行了离散傅里叶变换后的点值 { d i } \{d_i\} {di​}，即

d k   =   ∑ i = 0 n − 1 a i × ω n i k d_k~=~\sum_{i = 0}^{n - 1} a_i \times \omega_n^{ik} dk​ = i=0∑n−1​ai​×ωnik​

这个过程为DFT，使用FFT在 O ( n l o g n ) O(nlogn) O(nlogn)的时间复杂度下完成这一操作。

现在试图还原系数数列 { a i } \{a_i\} {ai​}。

结论：

a k   =   1 n ∑ i = 0 n − 1 d i ω n − k i a_k~=~\frac{1}{n} \sum_{i = 0}^{n - 1} d_i \omega_n^{-ki} ak​ = n1​i=0∑n−1​di​ωn−ki​

这个过程为IDFT，同样可以FFT在 O ( n l o g n ) O(nlogn) O(nlogn)的时间复杂度下完成这一操作。

如何使用FFT来计算 B ( x )   =   ∑ i = 0 n − 1 b i × x i B B(x)~=~\sum_{i = 0}^{n - 1} b_i \times x^iB B(x) = ∑i=0n−1​bi​×xiB

在 w n − k i w_n^{-ki} wn−ki​ ，其中 0 ≤ k < n 0 \leq k < n 0≤k<n 处的点值。

而 w n − k w_n^{-k} wn−k​ 可以看做 n n n 次本原单位根每次逆时针旋转本原单位根幅角的弧度，因此 ω n − k \omega_n^{-k} ωn−k​ 和 ω n k \omega_n^k ωnk​是一一对应的。具体的， w n − k = w n k + n w_n^{-k} = w_n^{k + n} wn−k​=wnk+n​ 。因此我们只需要使用 FFT 的方法，求出 B ( x ) B B(x)B B(x)B 在 ω n \omega_n ωn​ 各个幂次下的值，然后数组反过来，即令 a k   =   1 n ∑ i = 0 n B ( w n n − k ) a_k~=~\frac{1}{n} \sum_{i = 0}^n B(w_n^{n - k}) ak​ = n1​∑i=0n​B(wnn−k​) 即可。

这一步快速计算插值的过程叫做快速傅里叶逆变换 ( I n v e r s e   F a s t   F o u r i e r   T r a n s f o r m   ,   I F F T ) \tt (Inverse\ Fast\ Fourier\ Transform\ ,\ IFFT) (Inverse Fast Fourier Transform , IFFT)。

至此，我们得到了一个时间复杂度为 O ( n log ⁡ n ) O(n \log n) O(nlogn) 的多项式乘法计算方法。

算法流程：

在这里插入图片描述

0x21.2. 递归版

我们把上述的思路实现用递归可以很形象地实现FFT函数。

这里要注意一下，尽管C++自带有复数库，但是建议手写一下，代码不长，手写常数小。

void FFT(int limit, Complex *a, int type) {
    if (limit == 1) return ; //只有一个常数项
    Complex a1[limit >> 1], a2[limit >> 1];
    for (int i = 0; i <= limit; i += 2) //根据下标的奇偶性分类
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
    FFT(limit >> 1, a1, type);
    FFT(limit >> 1, a2, type);
    Complex Wn = Complex(cos(2.0 * Pi / limit) , type * sin(2.0 * Pi / limit)), w = Complex(1, 0);
    //Wn为单位根，w表示幂
    for (int i = 0; i < (limit >> 1); i++, w = w * Wn) //这里的w相当于公式中的k
        a[i] = a1[i] + w * a2[i],
               a[i + (limit >> 1)] = a1[i] - w * a2[i]; //利用单位根的性质，O(1)得到另一部分
}

int main() {
    int N = read(), M = read();
    for (int i = 0; i <= N; i++) a[i].x = read();
    for (int i = 0; i <= M; i++) b[i].x = read();
    int limit = 1; while (limit <= N + M) limit <<= 1;
    FFT(limit, a, 1);
    FFT(limit, b, 1);
    //后面的1表示要进行的变换是什么类型
    //1表示从系数变为点值
    //-1表示从点值变为系数
    //至于为什么这样是对的，可以参考一下c向量的推导过程，
    for (int i = 0; i <= limit; i++)
        a[i] = a[i] * b[i];
    FFT(limit, a, -1);
    for (int i = 0; i <= N + M; i++) printf("%d ", (int)(a[i].x / limit + 0.5)); //按照我们推倒的公式，这里还要除以n
    return 0;
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31

但是这样写常数太大了，还需要一个很大的数组存，实测会T飞

0x21.3. 迭代版（蝴蝶变换）

所以我们引入一种迭代版，利用蝴蝶变换巧妙地实现FFT，使其真正成为 O ( n l o g n ) O(nlogn) O(nlogn)的高效算法。

P3803 【模板】多项式乘法（FFT）

在这里插入图片描述

const double PI = acos(-1);
int n, m;
int res, ans[N];
int limit = 1;//补齐的2的整数幂N
int L;//二进制的位数
int R[N];//二进制翻转
struct Complex
{
    double x, y;
    Complex (double x = 0, double y = 0) : x(x), y(y) { }
}a[N], b[N];
//复数乘法：模长相乘，幅度相加
Complex operator * (Complex J, Complex Q) {return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);}
Complex operator - (Complex J, Complex Q) {return Complex(J.x - Q.x, J.y - Q.y);}
Complex operator + (Complex J, Complex Q) {return Complex(J.x + Q.x, J.y + Q.y);}

void FFT(Complex * A, int type)
{
    for(int i = 0; i < limit; ++ i)
        if(i < R[i])
            swap(A[i], A[R[i]]);
        //i小于R[i]时才交换，防止同一个元素交换两次，回到它原来的位置。
    
    //从底层往上合并
    for(int mid = 1; mid < limit; mid <<= 1) {
        //待合并区间长度的一半，最开始是两个长度为1的序列合并,mid = 1;
        Complex wn(cos(PI / mid), type * sin(PI / mid));//单位根w_n^1;

        for(int len = mid << 1, pos = 0; pos < limit; pos += len) {
            //len是区间的长度，pos是当前的位置,也就是合并到了哪一位
            Complex w(1, 0);//幂,一直乘，得到平方，三次方...

            for(int k = 0; k < mid; ++ k, w = w * wn) {
                //只扫左半部分，蝴蝶变换得到右半部分的答案,w 为 w_n^k
                Complex x = A[pos + k];//左半部分
                Complex y = w * A[pos + mid + k];//右半部分
                A[pos + k] = x + y;//左边加
                A[pos + mid + k] = x - y;//右边减
            }
        }
    }
    if(type == 1) return ;
    for(int i = 0; i <= limit; ++ i)
        A[i].x /= limit;
        //最后要除以limit也就是补成了2的整数幂的那个N，将点值转换为系数
        //（前面推过了点值与系数之间相除是N）
}
int main()
{
    n = read(), m = read();
    //读入多项式的每一项，保存在复数的实部
    for(int i = 0; i <= n; ++ i)
        a[i].x = read();
    for(int i = 0; i <= m; ++ i)
        b[i].x = read();
    while(limit <= n + m)
        limit <<= 1, L ++ ;
    //也可以写成：limit = 1 << int(log2(n + m) + 1); 
    // 补成2的整次幂，也就是N
    for(int i = 0; i < limit; ++ i)
        R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
    FFT(a, 1);//FFT 把a的系数表示转化为点值表示
    FFT(b, 1);//FFT 把b的系数表示转化为点值表示
    //计算两个系数表示法的多项式相乘后的点值表示
    for(int i = 0; i <= limit; ++ i)
        a[i] = a[i] * b[i];
    //对应项相乘，O(n)得到点值表示的多项式的解C，利用逆变换完成插值得到答案C的点值表示
    FFT(a, -1);

    for(int i = 0; i <= n + m; ++ i)
    	//这里的 x 和 y 是 double 的 hhh
        printf("%d ", (int)(a[i].x + 0.5));//注意要+0.5，否则精度会有问题
}


  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74

0x21.4. “三步变两步”优化

设 P P P和 Q Q Q 是 实多项式， F = P + Q i F = P+Qi F=P+Qi ，则 F 2 = P 2 − Q 2 + 2 P Q i F^2 = P^2-Q^2+2PQi F2=P2−Q2+2PQi ，注意到我们要求的 P Q PQ PQ正是 F F F 虚部的一半。这样只需要两次FFT就可以求出结果。

void FFT(Complex * A, int type)//FFT板子
{
    for(int i = 0; i < limit; ++ i)
        if(i < R[i])
            swap(A[i], A[R[i]]);

    for(int mid = 1; mid < limit; mid <<= 1) {
        Complex wn(cos(PI / mid), type * sin(PI / mid));

        for(int len = mid << 1, pos = 0; pos < limit; pos += len) {
            Complex w(1, 0);

            for(int k = 0; k < mid; ++ k, w = w * wn) {
                Complex x = A[pos + k];
                Complex y = w * A[pos + mid + k];
                A[pos + k] = x + y;
                A[pos + mid + k] = x - y;
            }
        }
    }
    if(type == 1) return ;
    for(int i = 0; i <= limit; ++ i)
        A[i].x /= limit, a[i].y /= limit;
}

int main()
{
    n = read(), m = read();
    for(int i = 0; i <= n; ++ i)
        a[i].x = read();
    for(int i = 0; i <= m; ++ i)
        a[i].y = read();//把b(x)放到a(x)的虚部上
    while(limit <= n + m)
        limit <<= 1, L ++ ;
    for(int i = 0; i < limit; ++ i)
        R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
    FFT(a, 1);
    for(int i = 0; i <= limit; ++ i)
        a[i] = a[i] * a[i];//求出a(x)^2
    FFT(a, -1);

    for(int i = 0; i <= n + m; ++ i)
        printf("%d ", (int)(a[i].y / 2 + 0.5)); 
        //虚部取出来除2，注意要+0.5，否则精度会有问题,这里的x和y都是double

}


  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47

0x22.快速数论变换（NTT）

0x22.1NTT简述

在这里插入图片描述

在这里插入图片描述

注：这里的 a − i k a^{-ik} a−ik指的是 a i k a^{ik} aik在模M意义下的逆元

模板题目大意：输入 n n n 次多项式 f ( x ) f(x) f(x) 和 m m m 次多项式 g ( x ) g(x) g(x)，求两个多项式的卷积

这里我们把多项式乘法封装成函数，更容易调用（NTT是多项式问题中最常用的变换）

const int p = 998244353, G = 3, Gi = 332748118;//这里的Gi是G的除法逆元
const int N = 5000007;

const double PI = acos(-1);

int n, m;
int res, ans[N];
int limit = 1;//
int L;//二进制的位数
int RR[N];
ll a[N], b[N];

inline int read()
{
    register int x = 0, f = 1;
    register char ch = getchar();
    while(ch < '0' || ch > '9') {if(ch == '-')f = -1;ch = getchar();}
    while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0';ch = getchar();}
    return x * f;
}

ll qpow(ll a, ll b)
{
    ll res = 1;
    while(b) {
        if(b & 1) res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res % p;
}

ll inv(ll x) {return qpow(x, p - 2);}

void NTT(ll *A, int type)
{
    for(int i = 0; i < limit; ++ i)
        if(i < RR[i])
            swap(A[i], A[RR[i]]);
    for(int mid = 1; mid < limit; mid <<= 1) {//原根代替单位根
        //ll wn = qpow(type == 1 ? G : Gi, (p - 1) / (mid << 1));
        ll wn = qpow(G, (p - 1) / (mid * 2));
        if(type == -1) wn = qpow(wn, p - 2);
        //如果超时了上面if这句话删掉，在下面的if(type == -1)里加上下面这个循环
        /*for (int i = 1; i < limit / 2; i ++)
        swap(A[i], A[limit - i]); */
        //逆变换则乘上逆元,因为我们算出来的公式中逆变换是(a^-ij)，也就是(a^ij)的逆元
        for(int len = mid << 1, pos = 0; pos < limit; pos += len) {
            ll w = 1;
            for(int k = 0; k < mid; ++ k, w = (w * wn) % p) {
                int x = A[pos + k], y = w * A[pos + mid + k] % p;
                A[pos + k] = (x + y) % p;
                A[pos + k + mid] = (x - y + p) % p;

            }
        }
    }

    if(type == -1) {
        ll limit_inv = inv(limit);//N的逆元（N是limit, 指的是2的整数幂）
        for(int i = 0; i < limit; ++ i)
            A[i] = (A[i] * limit_inv) % p;//NTT还是要除以n的，但是这里把除换成逆元了，inv就是n在模p意义下的逆元
    }
}//代码实现上和FFT相差无几
//多项式乘法
void poly_mul(ll *a, ll *b, int deg)
{
    for(limit = 1, L = 0; limit <= deg; limit <<= 1) L ++ ;
    for(int i = 0; i < limit; ++ i) {
        RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
    }
    NTT(a, 1);
    NTT(b, 1);
    for(int i = 0; i < limit; ++ i) a[i] = a[i] * b[i] % p;
    NTT(a, -1);
}

int main()
{
    n = read(), m = read();
    for(int i = 0; i <= n; ++ i) a[i] = (read() + p) % p;//取模好习惯
    for(int i = 0; i <= m; ++ i) b[i] = (read() + p) % p;
    poly_mul(a, b, n + m);
    for(int i = 0; i <= n + m; ++ i)
        printf("%d ", a[i]);
    return 0;
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87

下面的所有内容待更（待学），预计2021年前更完hhh

0x23.快速沃尔什变换（FWT）

《小学生都能看懂的快速沃尔什变换从入门到升天教程》（FWT / FMT / FMI）（最最严谨清晰的证明！零基础也能得学会！）

0x24.快速莫比乌斯变换（FMT）

0x25.任意模数变换（MTT）

P4239 任意模数多项式乘法逆

算法学习FFT系列（4）：任意模数的快速傅里叶变换（MTT）

0x30.多项式求逆

0x31. 多项式求逆简述

0x32. 多项式求逆模板

P4238 【模板】多项式乘法逆

在这里插入图片描述
直接按照上述的思路实现一下模板即可。注意要中间 h ( x ) h(x) h(x) 要取的是模 x ⌈ n ⌉ 2 x^{\frac{\lceil n \rceil}{2}} x2⌈n⌉​ 意义下的 g ( x ) g(x) g(x) ，所以我们应该把后面大于 d e g deg deg 的全部置为 0 0 0 （这里的 d e g deg deg 实际上就是我们证明的时候用到的 n n n ），我们可以直接使用 f i l l fill fill 函数实现这一操作。

const int N = 5000007;
const int p = 998244353, g = 3, gg = 3, ig = 332738118;
int mod = p;
int n, m;
int limit = 1;
int L;
int R[N];
ll A[N], B[N];

template <typename T>void read(T &x)
{
    x = 0;
    register int f = 1;
    register char ch = getchar();
    while(ch < '0' || ch > '9') {if(ch == '-')f = -1;ch = getchar();}
    while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0';ch = getchar();}
    x *= f;
}

ll qpow(ll a, ll b)
{
    ll res = 1;
    while(b) {
        if(b & 1) res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res % p;
}

ll inv(ll x) {return qpow(x, p - 2);}

void NTT(ll *A, int type)
{
    for(int i = 0; i < limit; ++ i) {
        if(i < R[i])
            swap(A[i], A[R[i]]);
    }
    for(int mid = 1; mid < limit; mid <<= 1) {
        ll wn = qpow(g, (p - 1) / (mid * 2));
        if(type == -1) wn = qpow(wn, p - 2);
		//如果超时了上面if这句话删掉，在下面的if(type == -1)里加上下面这个循环
        /*for (int i = 1; i < limit / 2; i ++)
        swap(A[i], A[limit - i]); */
        for(int len = mid << 1, pos = 0; pos < limit; pos += len) {
            ll w = 1;
            for(int k = 0; k < mid; ++ k, w = (w * wn) % p) {
                int x = A[pos + k] % p, y = w * A[pos + mid + k] % p;
                A[pos + k] = (x + y) % p;
                A[pos + k + mid] = (x - y + p) % p;
            }
        }
    }
    if(type == 1) return ;
    ll inv_limit = inv(limit);
    for(int i = 0; i < limit; ++ i)
        A[i] = (A[i] * inv_limit) % p;
}

//。C[N] == f[N] % x ^ deg;
//。B[N] == h[N] == g[N] % x ^ ⌈deg/2⌉
//。A[N] == f[N];
ll C[N];

void get_inv(ll *A, ll *B, int deg)//deg = 多项式的度
{
    if(deg == 1) {
        B[0] = inv(A[0]);//A[0]的逆即常数项
        return ;
    }

    get_inv(A, B, (deg + 1) >> 1);//递归分治
    for(limit = 1; limit <= (deg << 1); limit <<= 1);
    for(int i = 0; i < limit; ++ i) {
        R[i] = (R[i >> 1] >> 1) | ((i & 1) ? (limit >> 1) : 0);
        C[i] = (i < deg ? A[i] : 0);//只算到⌈n/2⌉，后面的全部为0;
    }
    NTT(C, 1), NTT(B, 1);
    for(int i = 0; i < limit; ++ i) {
        B[i] = (2ll - C[i] * B[i] % p + p) % p * B[i] % p;
    }
    NTT(B, -1);
    fill(B + deg, B + limit, 0);//非常重要，因为是在模 x^deg 意义下 ，所以大于deg的置0
}

int main()
{
    read(n);
    for(int i = 0; i < n; ++ i) read(A[i]);
    get_inv(A, B, n);
    for(int i = 0; i < n; ++ i) {
        printf("%lld ", B[i]);
    }
    return 0;
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91
  • 92
  • 93
  • 94
  • 95

0x40.多项式开方

0x41. 多项式开方简述

在这里插入图片描述

0x42. 多项式开方(根)模板

在这里插入图片描述
在这里插入图片描述

const int N = 5000007;
const int p = 998244353, gg = 3, ig = 332738118;
const int mod = 998244353;
int limit = 1;
int L;
int R[N];
ll f[N], g[N], h[N];
ll A[N], B[N], C[N];
//部分模板代码省略
template <typename T>void read(T &x);//快读
ll qpow(ll a, ll b);//快速幂
ll inv(ll x) {return qpow(x, mod - 2);}//逆元
int inv2 = qpow(2, mod - 2);
void NTT(ll *A, int type);//快速数论变换
void get_inv(ll *A, ll *B, int deg);//多项式求逆
void poly_sqrt(ll *f, ll *h, int deg) {//polynomial 多项式
    if (deg == 1) {
        h[0] = 1;
        return;
    }

    poly_sqrt(f, h, deg + 1 >> 1);

    limit = 1;
    while (limit < deg << 1) {
        limit <<= 1;
    }
    fill(g, g + limit, 0);
    get_inv(h, g, deg);
    copy(f, f + deg, t);
    fill(t + deg, t + limit, 0);
    NTT(t, 1);
    NTT(g, 1);
    NTT(h, 1);

    for (int i = 0; i < limit; i++) {
        h[i] = 1LL * inv2 * (1LL * h[i] % mod + 1LL * g[i] * t[i] % mod) % mod;
    }
    NTT(h, -1);
    fill(h + deg, h + limit, 0);
}
int n, m;
ll a[N], b[N];

//输入g输出f中间多项式h
//h是g的一半的平方根
int main()
{
    read(n);
    for(int i = 0; i < n; ++ i) read(a[i]);

    poly_sqrt(a, b, n);

    for(int i = 0; i < n; ++ i)
        printf("%lld ", b[i]);
    return 0;
}



  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59

0x50.多项式除法 / 取模

0x51.多项式除法简述

题目大意：给定一个 n n n次多项式 F ( x ) F(x) F(x)， m m m次多项式 G ( x ) G(x) G(x)，求一个 n − m n-m n−m次多项式 Q ( x ) Q(x) Q(x)，和一个小于 m m m次的多项式 R ( x ) R(x) R(x)。

满足 F ( x ) = G ( x ) Q ( x ) + R ( x ) F(x)=G(x)Q(x)+R(x) F(x)=G(x)Q(x)+R(x)

即求 F ( x ) F(x) F(x)除以 G ( x ) G(x) G(x)得到的商 Q ( x ) Q(x) Q(x)以及余数 R ( x ) R(x) R(x)。

定义 F r ( x ) Fr(x) Fr(x)为多项式 F ( x ) F(x) F(x)系数翻转后的多项式，即

F ( x ) = ∑ i = 0 n a i x i F(x)=\sum^n_{i=0}a_ix^i F(x)=∑i=0n​ai​xi， G ( x ) = ∑ i = 0 n a n − i x i G(x)=\sum^n_{i=0}a_{n-i}x^i G(x)=∑i=0n​an−i​xi

翻转的实现可以使用公式 F r ( x ) = x n F ( x − 1 ) → F ( x ) = x n F r ( x − 1 ) Fr(x)=x^nF(x^{-1}) \to F(x)=x^nFr(x^{-1}) Fr(x)=xnF(x−1)→F(x)=xnFr(x−1)，也可以直接循环来交换系数来实现

最终推导得到： F r ( x ) ≡ Q ( x ) G r ( x )   ( m o d   x n − m + 1 ) Fr(x)\equiv Q(x)Gr(x)\ (mod\ x^{n-m+1}) Fr(x)≡Q(x)Gr(x) (mod xn−m+1)

即：

Q r ( x ) ≡ F r ( x ) G r − 1 ( x )   ( m o d   x n − m + 1 ) Qr(x)\equiv Fr(x)Gr^{-1}(x)\ (mod\ x^{n-m+1}) Qr(x)≡Fr(x)Gr−1(x) (mod xn−m+1)

R ( x ) = F ( x ) − G ( x ) Q ( x ) R(x)=F(x)-G(x)Q(x) R(x)=F(x)−G(x)Q(x)

简单实现一下即可。

时间复杂度： O ( n l o g n ) O(nlogn) O(nlogn)

0x52. 多项式除法模板

const int N = 4000007;
const int p = 998244353, gg = 3, ig = 332738118;
const int mod = 998244353;
int limit = 1;
int L;
int RR[N];
ll F[N], G[N], H[N], Q[N], R[N];
ll A[N], B[N], C[N];
ll gg_inv;

template <typename T>void read(T &x)
{
    x = 0;
    register int f = 1;
    register char ch = getchar();
    while(ch < '0' || ch > '9') {if(ch == '-')f = -1;ch = getchar();}
    while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0';ch = getchar();}
    x *= f;
}

ll qpow(ll a, ll b)
{
    ll res = 1;
    while(b) {
        if(b & 1) res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res % p;
}

ll inv(ll x) {return qpow(x, mod - 2);}

void NTT(ll *A, int type)
{
    for(int i = 0; i < limit; ++ i)
        if(i < RR[i])
            swap(A[i], A[RR[i]]);
    for(int mid = 1; mid < limit; mid <<= 1) {
        ll wn = qpow(gg, (mod - 1) / (2 * mid));
        if(type == -1) wn = qpow(wn, mod - 2);
        for(int len = mid << 1, pos = 0; pos < limit; pos += len) {
            ll w = 1;
            for(int k = 0; k < mid; ++ k, w = (w * wn) % mod){
                int x = A[pos + k] % mod;
                int y = w * A[pos + k + mid] % mod;
                A[pos + k] = (x + y) % mod;
                A[pos + k + mid] = (x - y + mod) % mod;
            }
        }
    }
    if(type == -1) {
        ll limit_inv = inv(limit);
        for(int i = 0; i < limit; ++ i) {
            A[i] = (A[i] * limit_inv) % mod;
        }
    }
}

//多项式求逆
void get_inv(ll *A, ll *B, int deg)
{
    if(deg == 1) {
        B[0] = inv(A[0]);
        return ;
    }

    get_inv(A, B, (deg + 1) >> 1);

    for(limit = 1, L = 0; limit <= (deg << 1) * 2; limit <<= 1) L ++ ;
    for(int i = 0; i < limit; ++ i) {
        RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
        C[i] = (i < deg ? A[i] : 0);
    }
    NTT(C, 1), NTT(B, 1);
    for(int i = 0; i <= limit; ++ i)
        B[i] = (2ll - C[i] * B[i] % mod + mod) % mod * B[i] % mod;
    NTT(B, -1);
    fill(B + deg, B + limit, 0);
}

void poly_mul(ll *a, ll *b, int n, int m)
{
    for(limit = 1, L = 0; limit <= (n + m) * 2; limit <<= 1) L ++ ;
    for(int i = 0; i < limit; ++ i) {
        RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
    }
    NTT(a, 1);
    NTT(b, 1);
    for(int i = 0; i < limit; ++ i ) {
        a[i] = a[i] * b[i] % mod;
    }
    NTT(a, -1);
}

ll t[N], Fr[N], Gr[N], Qr[N], Gr_inv[N];
int n, m;

int main()
{
    //gg_inv = inv(gg);
    read(n);
    read(m);
    //n ++ , m ++ ;
    for(int i = 0; i <= n; ++ i) read(F[i]), Fr[n - i] = F[i];//读入F，求翻转的Fr数组（可以直接手动翻转hhh）
    for(int i = 0; i <= m; ++ i) read(G[i]), Gr[m - i] = G[i];//读入G，求翻转的Gr数组
    for(int i = n - m + 2; i <= m; ++ i) Gr[i] = 0;//Gr数组应该只有n - m项
    get_inv(Gr, Gr_inv, n - m + 1);//求逆

    poly_mul(Fr, Gr_inv, n, n - m);//乘起来得到Qr

    for(int i = 0; i <= n - m; ++ i) Q[i] = Fr[n - m - i];//求Qr数组的翻转既是要求的Q数组（商）

    for(int i = 0; i <= n - m; ++ i) printf("%lld ", Q[i]);

    puts("");

    poly_mul(G, Q, m, n - m);
    //R(x) = F(x) - G(x)Q(x);小于m项
    for(int i = 0; i < m; ++ i) R[i] = (F[i] - G[i] + mod) % mod;

    for(int i = 0; i < m; ++ i)
        printf("%lld ", R[i]);

    puts("");

    return 0;
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91
  • 92
  • 93
  • 94
  • 95
  • 96
  • 97
  • 98
  • 99
  • 100
  • 101
  • 102
  • 103
  • 104
  • 105
  • 106
  • 107
  • 108
  • 109
  • 110
  • 111
  • 112
  • 113
  • 114
  • 115
  • 116
  • 117
  • 118
  • 119
  • 120
  • 121
  • 122
  • 123
  • 124
  • 125
  • 126
  • 127
  • 128

vector 版本还没有调出来呜呜呜

0x60.多项式牛顿迭代

0x61. 多项式牛顿迭代及其应用

在这里插入图片描述

0x70.多项式多点求值 / 快速插值

0x80.多项式三角函数

【学习笔记】超简单的多项式三角函数（含全套证明）

g ( x ) ≡ cos ⁡ f ( x ) ≡ e i ⋅ f ( x ) + e − i ⋅ f ( x ) 2 (   m o d     x n ) g(x) \equiv \cos f(x) \equiv \frac{e^{i \cdot f(x)}+e^{-i \cdot f(x)}}{2}\left(\bmod\ x^{n}\right) g(x)≡cosf(x)≡2ei⋅f(x)+e−i⋅f(x)​(mod xn)

g ( x ) ≡ sin ⁡ f ( x ) ≡ e i ⋅ f ( x ) − e − i ⋅ f ( x ) 2 i (   m o d     x n ) g(x) \equiv \sin f(x) \equiv \frac{e^{i \cdot f(x)}-e^{-i \cdot f(x)}}{2 i}\left(\bmod\ x^{n}\right) g(x)≡sinf(x)≡2iei⋅f(x)−e−i⋅f(x)​(mod xn)

其中：

i 2 ≡ − 1 ≡ 998244352 ( m o d 998244353 ) i^{2} \equiv-1 \equiv 998244352(\mathrm{mod} 998244353) i2≡−1≡998244352(mod998244353)

较小解为：

i = 86583718 i=86583718 i=86583718

0x90.多项式反三角函数

【学习笔记】超简单的多项式反三角函数（含全套证明）

反正弦函数

g ( x ) ≡ ∫ f ′ ( x ) 1 − f 2 ( x ) d x (   m o d     x n ) g(x) \equiv \int \frac{f^{\prime}(x)}{\sqrt{1-f^{2}(x)}} d x \quad\left(\bmod\ x^{n}\right) g(x)≡∫1−f2(x) ​f′(x)​dx(mod xn)
反正切函数

g ( x ) ≡ ∫ f ′ ( x ) 1 + f 2 ( x ) d x (   m o d     x n ) g(x) \equiv \int \frac{f^{\prime}(x)}{1+f^{2}(x)} d x \quad\left(\bmod\ x^{n}\right) g(x)≡∫1+f2(x)f′(x)​dx(mod xn)

0x100.常系数齐次线性递推（多项式优化）

0x110.快速差分/前缀和

0x120.多项式对数（多项式求导 + 积分）

0x121. 多项式对数简述（含常见函数求导表）

其实非常简单，就是直接求导再积分即可。
在这里插入图片描述

0x122. 多项式对数函数（多项式 ln）模板

在这里插入图片描述

直接简单实现一下就行了。

我们输入多项式 f f f ，先求 f f f的导数 A A A ， f f f 的逆元 B B B，然后把他们乘起来（ N T T NTT NTT），最后求积分 g g g 既是答案，输出即可。

至于求导和积分，因为我们的多项式 f ( x ) = a 0 + a 1 × x + a 2 × x 2 … a n × x n − 1 f(x) = a_0 + a_1 \times x + a_2 \times x^2 \dots a_n \times x^{n - 1} f(x)=a0​+a1​×x+a2​×x2…an​×xn−1

( x a ) ′ = a x a − 1   ,   ∫ x a d x = 1 a + 1 x a + 1 (x^{a})'=ax^{a-1}\ ,\ \int x^adx=\frac{1}{a+1}x^{a+1} (xa)′=axa−1 , ∫xadx=a+11​xa+1

const int N = 5000007;
const int p = 998244353, gg = 3, ig = 332738118;
const int mod = 998244353;
int limit = 1;
int L;
int RR[N];
ll f[N], g[N];
ll A[N], B[N], C[N];
//省略部分上面已有的代码
template <typename T>void read(T &x);//快读
ll qpow(ll a, ll b);//快速幂
ll inv(ll x) {return qpow(x, mod - 2);}//求乘法逆元
void NTT(ll *A, int type);//快速数论变换
void get_inv(ll *A, ll *B, int deg);//多项式求逆
//多项式求导
void get_dev(ll *A, ll *B, int n)
{
    for(int i = 1; i < n; ++ i)
        B[i - 1] = i * A[i] % mod;
    B[n - 1] = 0;
}
//多项式求积分
void get_indev(ll *A, ll *B, int n)
{
    for(int i = 1; i < n; ++ i)
        B[i] = A[i - 1] * qpow(i, mod - 2) % mod;
    B[0] = 0;
}

//多项式乘法
void poly_mul(ll *a, ll *b, int deg)
{
    for(limit = 1, L = 0; limit <= deg; limit <<= 1) L ++ ;
    for(int i = 0; i < limit; ++ i) {
        RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
    }
    NTT(a, 1);
    NTT(b, 1);
    for(int i = 0; i < limit; ++ i) a[i] = a[i] * b[i] % mod;
    NTT(a, -1);
}
//多项式对数
void get_ln(ll *f, ll *g, int n)
{
    get_dev(f, A, n);
    get_inv(f, B, n);
    poly_mul(A, B, n);
    get_indev(A, g, n);
}

int n, m;

int main()
{
    read(n);
    for(int i = 0; i < n; ++ i) read(f[i]);
    for(limit = 1; limit <= n; limit <<= 1);
    get_ln(f, g, limit);

    for(int i = 0; i < n; ++ i)
        printf("%d ", g[i]);
    puts("");
    return 0;
}


  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65

0x130.多项式指数（大毒瘤 ）

由牛顿迭代得到的指数公式：

f ( x ) ≡ f ( x 0 ) ( 1 − l n ( f ( x 0 ) ) + g ( x ) )   m o d (   x n ) f(x) \equiv f(x_0)(1-ln(f(x_0))+g(x))\ mod(\ x^n) f(x)≡f(x0​)(1−ln(f(x0​))+g(x)) mod( xn)

代码见下面的全家桶

0x140.多项式快速幂

我们知道指数有一个性质：

l n x n = n   l n x lnx^n=n\ lnx lnxn=n lnx

我们要求的是 f ( x ) k f(x)^k f(x)k

可以很自然地构造出：

f ( x ) k = e l n f ( x ) k = e k   l n f ( x ) f(x)^k=e^{lnf(x)^k}=e^{k\ lnf(x)} f(x)k=elnf(x)k=ek lnf(x)

所以我们只需要求 f ( x ) f(x) f(x)的 对数ln 乘上 k k k然后求 指数 就行了QwQ

代码见下面的全家桶

0x150.多项式全家桶！！！（code）

“优 雅”

//#pragma GCC optimize(2)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N = 3000007;
const int p = 998244353, gg = 3, ig = 332738118, img = 86583718;
const int mod = 998244353;

template <typename T>void read(T &x)
{
    x = 0;
    register int f = 1;
    register char ch = getchar();
    while(ch < '0' || ch > '9') {if(ch == '-')f = -1;ch = getchar();}
    while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0';ch = getchar();}
    x *= f;
}

int qpow(int a, int b)
{
    int res = 1;
    while(b) {
        if(b & 1) res = 1ll * res * a % mod;
        a = 1ll * a * a % mod;
        b >>= 1;
    }
    return res;
}

namespace Poly
{
    #define mul(x, y) (1ll * x * y >= mod ? 1ll * x * y % mod : 1ll * x * y)
    #define minus(x, y) (1ll * x - y < 0 ? 1ll * x - y + mod : 1ll * x - y)
    #define plus(x, y) (1ll * x + y >= mod ? 1ll * x + y - mod : 1ll * x + y)
    #define ck(x) (x >= mod ? x - mod : x)//取模运算太慢了

    typedef vector<int> poly;
    const int G = 3;//根据具体的模数而定，原根可不一定不一样！！！
    //一般模数的原根为 2 3 5 7 10 6
    const int inv_G = qpow(G, mod - 2);
    int RR[N], deer[2][19][N], inv[N];

    void init(const int t) {//预处理出来NTT里需要的w和wn，砍掉了一个log的时间
        for(int p = 1; p <= t; ++ p) {
            int buf1 = qpow(G, (mod - 1) / (1 << p));
            int buf0 = qpow(inv_G, (mod - 1) / (1 << p));
            deer[0][p][0] = deer[1][p][0] = 1;
            for(int i = 1; i < (1 << p); ++ i) {
                deer[0][p][i] = 1ll * deer[0][p][i - 1] * buf0 % mod;//逆
                deer[1][p][i] = 1ll * deer[1][p][i - 1] * buf1 % mod;
            }
        }
        inv[1] = 1;
        for(int i = 2; i <= (1 << t); ++ i)
            inv[i] = 1ll * inv[mod % i] * (mod - mod / i) % mod;
    }

    int NTT_init(int n) {//快速数论变换预处理
        int limit = 1, L = 0;
        while(limit <= n) limit <<= 1, L ++ ;
        for(int i = 0; i < limit; ++ i)
            RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
        return limit;
    }

    void NTT(poly &A, int type, int limit) {//快速数论变换
        A.resize(limit);
        for(int i = 0; i < limit; ++ i)
            if(i < RR[i])
                swap(A[i], A[RR[i]]);
        for(int mid = 2, j = 1; mid <= limit; mid <<= 1, ++ j) {
            int len = mid >> 1;
            for(int pos = 0; pos < limit; pos += mid) {
                int *wn = deer[type][j];
                for(int i = pos; i < pos + len; ++ i, ++ wn) {
                    int tmp = 1ll * (*wn) * A[i + len] % mod;
                    A[i + len] = ck(A[i] - tmp + mod);
                    A[i] = ck(A[i] + tmp);
                }
            }
        }
        if(type == 0) {
            for(int i = 0; i < limit; ++ i)
                A[i] = 1ll * A[i] * inv[limit] % mod;
        }
    }

    poly poly_mul(poly A, poly B) {//多项式乘法
        int deg = A.size() + B.size() - 1;
        int limit = NTT_init(deg);
        poly C(limit);
        NTT(A, 1, limit);
        NTT(B, 1, limit);
        for(int i = 0; i < limit; ++ i)
            C[i] = 1ll * A[i] * B[i] % mod;
        NTT(C, 0, limit);
        C.resize(deg);
        return C;
    }

    poly poly_inv(poly &f, int deg) {//多项式求逆
        if(deg == 1)
            return poly(1, qpow(f[0], mod - 2));

        poly A(f.begin(), f.begin() + deg);
        poly B = poly_inv(f, (deg + 1) >> 1);
        int limit = NTT_init(deg << 1);
        NTT(A, 1, limit), NTT(B, 1, limit);
        for(int i = 0; i < limit; ++ i)
            A[i] = B[i] * (2 - 1ll * A[i] * B[i] % mod + mod) % mod;
        NTT(A, 0, limit);
        A.resize(deg);
        return A;
    }

    poly poly_dev(poly f) {//多项式求导
        int n = f.size();
        for(int i = 1; i < n; ++ i) f[i - 1] = 1ll * f[i] * i % mod;
        return f.resize(n - 1), f;//f[0] = 0，这里直接扔了,从1开始
    }

    poly poly_idev(poly f) {//多项式求积分
        int n = f.size();
        for(int i = n - 1; i ; -- i) f[i] = 1ll * f[i - 1] * inv[i] % mod;
        return f[0] = 0, f;
    }

    poly poly_ln(poly f, int deg) {//多项式求对数
        poly A = poly_idev(poly_mul(poly_dev(f), poly_inv(f, deg)));
        return A.resize(deg), A;
    }

    poly poly_exp(poly &f, int deg) {//多项式求指数
        if(deg == 1)
            return poly(1, 1);

        poly B = poly_exp(f, (deg + 1) >> 1);
        B.resize(deg);
        poly lnB = poly_ln(B, deg);
        for(int i = 0; i < deg; ++ i)
            lnB[i] = ck(f[i] - lnB[i] + mod);

        int limit = NTT_init(deg << 1);//n -> n^2
        NTT(B, 1, limit), NTT(lnB, 1, limit);
        for(int i = 0; i < limit; ++ i)
            B[i] = 1ll * B[i] * (1 + lnB[i]) % mod;
        NTT(B, 0, limit);
        B.resize(deg);
        return B;
    }

    poly poly_sqrt(poly &f, int deg) {//多项式开方
        if(deg == 1) return poly(1, 1);
        poly A(f.begin(), f.begin() + deg);
        poly B = poly_sqrt(f, (deg + 1) >> 1);
        poly IB = poly_inv(B, deg);
        int limit = NTT_init(deg << 1);
        NTT(A, 1, limit), NTT(IB, 1, limit);
        for(int i = 0; i < limit; ++ i)
            A[i] = 1ll * A[i] * IB[i] % mod;
        NTT(A, 0, limit);
        for(int i =0; i < deg; ++ i)
            A[i] = 1ll * (A[i] + B[i]) * inv[2] % mod;
        A.resize(deg);
        return A;
    }

    poly poly_pow(poly f, int k) {//多项式快速幂
        f = poly_ln(f, f.size());
        for(auto &x : f) x = 1ll * x * k % mod;
        return poly_exp(f, f.size());
    }

    poly poly_cos(poly f, int deg) {//多项式三角函数（cos）
        poly A(f.begin(), f.begin() + deg);
        poly B(deg), C(deg);
        for(int i = 0; i < deg; ++ i)
            A[i] = 1ll * A[i] * img % mod;

        B = poly_exp(A, deg);
        C = poly_inv(B, deg);
        int inv2 = qpow(2, mod - 2);
        for(int i = 0; i < deg; ++ i)
            A[i] = 1ll * (1ll * B[i] + C[i]) % mod * inv2 % mod;
        return A;
    }

    poly poly_sin(poly f, int deg) {//多项式三角函数（sin）
        poly A(f.begin(), f.begin() + deg);
        poly B(deg), C(deg);
        for(int i = 0; i < deg; ++ i)
            A[i] = 1ll * A[i] * img % mod;

        B = poly_exp(A, deg);
        C = poly_inv(B, deg);
        int inv2i = qpow(img << 1, mod - 2);
        for(int i = 0; i < deg; ++ i)
            A[i] = 1ll * (1ll * B[i] - C[i] + mod) % mod * inv2i % mod;
        return A;
    }

    poly poly_arcsin(poly f, int deg) {
        poly A(f.size()), B(f.size()), C(f.size());
        A = poly_dev(f);
        B = poly_mul(f, f);
        for(int i = 0; i < deg; ++ i)
            B[i] = minus(mod, B[i]);
        B[0] = plus(B[0], 1);
        C = poly_sqrt(B, deg);
        C = poly_inv(C, deg);
        C = poly_mul(A, C);
        C = poly_idev(C);
        return C;
    }

    poly poly_arctan(poly f, int deg) {
        poly A(f.size()), B(f.size()), C(f.size());
        A = poly_dev(f);
        B = poly_mul(f, f);
        B[0] = plus(B[0], 1);
        C = poly_inv(B, deg);
        C = poly_mul(A, C);
        C = poly_idev(C);
        return C;
    }
}

using Poly::poly;
using Poly::poly_arcsin;
using Poly::poly_arctan;

int n, m, x, k, type;
poly f, g;
char s[N];

int main()
{
    Poly::init(18);//2^21 = 2,097,152,根据题目数据多项式项数的大小自由调整，注意大小需要跟deer数组同步(21+1=22)

    read(n), read(type);

    for(int i = 0; i < n; ++ i)
        read(x), f.push_back(x);

    if(type == 0) g = poly_arcsin(f, n);
    else g = poly_arctan(f, n);

    for(int i = 0; i < n; ++ i)
        printf("%d ", g[i]);
    return 0;
}

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91
  • 92
  • 93
  • 94
  • 95
  • 96
  • 97
  • 98
  • 99
  • 100
  • 101
  • 102
  • 103
  • 104
  • 105
  • 106
  • 107
  • 108
  • 109
  • 110
  • 111
  • 112
  • 113
  • 114
  • 115
  • 116
  • 117
  • 118
  • 119
  • 120
  • 121
  • 122
  • 123
  • 124
  • 125
  • 126
  • 127
  • 128
  • 129
  • 130
  • 131
  • 132
  • 133
  • 134
  • 135
  • 136
  • 137
  • 138
  • 139
  • 140
  • 141
  • 142
  • 143
  • 144
  • 145
  • 146
  • 147
  • 148
  • 149
  • 150
  • 151
  • 152
  • 153
  • 154
  • 155
  • 156
  • 157
  • 158
  • 159
  • 160
  • 161
  • 162
  • 163
  • 164
  • 165
  • 166
  • 167
  • 168
  • 169
  • 170
  • 171
  • 172
  • 173
  • 174
  • 175
  • 176
  • 177
  • 178
  • 179
  • 180
  • 181
  • 182
  • 183
  • 184
  • 185
  • 186
  • 187
  • 188
  • 189
  • 190
  • 191
  • 192
  • 193
  • 194
  • 195
  • 196
  • 197
  • 198
  • 199
  • 200
  • 201
  • 202
  • 203
  • 204
  • 205
  • 206
  • 207
  • 208
  • 209
  • 210
  • 211
  • 212
  • 213
  • 214
  • 215
  • 216
  • 217
  • 218
  • 219
  • 220
  • 221
  • 222
  • 223
  • 224
  • 225
  • 226
  • 227
  • 228
  • 229
  • 230
  • 231
  • 232
  • 233
  • 234
  • 235
  • 236
  • 237
  • 238
  • 239
  • 240
  • 241
  • 242
  • 243
  • 244
  • 245
  • 246
  • 247
  • 248
  • 249
  • 250
  • 251
  • 252

他的板子跑的实在是太快了

/*
namespace poly_base (变换基础)
namespace poly (多项式初等函数)
namespace CDQ_NTT (分治多项式技巧)
namespace poly_evaluation (多点求值)
namespace poly_interpolation (快速插值)
namespace linear_recur (常系数线性齐次递推式 - Fiduccia)
namespace linear_recur_single (常系数线性齐次递推式 - 单点)
namespace miscellaneous (杂项)
*/
// This is a NEW polynomial template for C++11
#include <bits/stdc++.h>
#define EB emplace_back
#define lg2 std::__lg

typedef unsigned long long u64;
const int N = 530000, mod = 998244353, iv2 = (mod + 1) / 2, unity = 31;
typedef int vec[N], *pvec;
typedef std::pair <int, int> pr;
typedef std::vector <int> vector;

vec inv, fact, finv;

inline int min(const int x, const int y) {return x < y ? x : y;}
inline int max(const int x, const int y) {return x < y ? y : x;}
inline int & reduce(int &x) {return x += x >> 31 & mod;}
inline int & neg(int &x) {return x = (!x - 1) & (mod - x);}
u64 PowerMod(u64 a, int n, u64 c = 1) {for (; n; n >>= 1, a = a * a % mod) if (n & 1) c = c * a % mod; return c;}

namespace poly_base {
	int l, n; u64 iv; vec w2;

	void init(int n = N, bool dont_calc_factorials = true) {
		int i, t;
		for (inv[1] = 1, i = 2; i < n; ++i) inv[i] = u64(mod - mod / i) * inv[mod % i] % mod;
		if (!dont_calc_factorials) for (*finv = *fact = i = 1; i < n; ++i) fact[i] = (u64)fact[i - 1] * i % mod, finv[i] = (u64)finv[i - 1] * inv[i] % mod;
		t = min(n > 1 ? lg2(n - 1) : 0, 21),
		*w2 = 1, w2[1 << t] = PowerMod(unity, 1 << (21 - t));
		for (i = t; i; --i) w2[1 << (i - 1)] = (u64)w2[1 << i] * w2[1 << i] % mod;
		for (i = 1; i < n; ++i) w2[i] = (u64)w2[i & (i - 1)] * w2[i & -i] % mod;
	}

	inline void NTT_init(int len) {n = 1 << (l = len), iv = mod - (mod - 1) / n;}

	void DIF(int *a) {
		int i, *j, *k, len = n >> 1, R, *o;
		for (i = 0; i < l; ++i, len >>= 1)
			for (j = a, o = w2; j != a + n; j += len << 1, ++o)
				for (k = j; k != j + len; ++k)
					R = (u64)*o * k[len] % mod, reduce(k[len] = *k - R), reduce(*k += R - mod);
	}

	void DIT(int *a) {
		int i, *j, *k, len = 1, R, *o;
		for (i = 0; i < l; ++i, len <<= 1)
			for (j = a, o = w2; j != a + n; j += len << 1, ++o)
				for (k = j; k != j + len; ++k)
					reduce(R = *k + k[len] - mod), k[len] = u64(*k - k[len] + mod) * *o % mod, *k = R;
	}

	inline void DNTT(int *a) {DIF(a);}
	inline void IDNTT(int *a) {
		DIT(a), std::reverse(a + 1, a + n);
		for (int i = 0; i < n; ++i) a[i] = a[i] * iv % mod;
	}

	inline void DIF(int *a, int *b) {memcpy(b, a, n << 2), DIF(b);}
	inline void DIT(int *a, int *b) {memcpy(b, a, n << 2), DIT(b);}
	inline void DNTT(int *a, int *b) {memcpy(b, a, n << 2), DNTT(b);}
	inline void IDNTT(int *a, int *b) {memcpy(b, a, n << 2), IDNTT(b);}
}

namespace poly {
	using namespace poly_base;

	vec B1, B2, B3, B4, B5, B6;

	// Multiplication (use one buffer, 3-dft of length 2n)
	void mul(int deg, pvec a, pvec b, pvec c) {
		if (!deg) {*c = (u64)*a * *b % mod; return;}
		NTT_init(lg2(deg) + 1), DNTT(a, c), DNTT(b, B1);
		for (int i = 0; i < n; ++i) c[i] = (u64)c[i] * B1[i] % mod;
		IDNTT(c);
	}

	// Inversion (use three buffers, 5-dft)
	void inv(int deg, pvec a, pvec b) {
		int i, len; assert(*a);
		if (*b = PowerMod(*a, mod - 2), deg <= 1) return;
		memset(b + 1, 0, i = 8 << lg2(deg - 1)), memset(B1, 0, i), *B1 = *a;

		for (len = 0; 1 << len < deg; ++len) {
			NTT_init(len + 1);

			memcpy(B1 + (n >> 1), a + (n >> 1), n << 1), DIF(b, B2), DIF(B1, B3);
			for (i = 0; i < n; ++i) B3[i] = (u64)B3[i] * B2[i] % mod; DIT(B3);
			for (i = n >> 1; i < n; ++i) B3[i] = B3[n - i] * iv % mod;

			memset(B3, 0, n << 1), DIF(B3);
			for (i = 0; i < n; ++i) B3[i] = (u64)B3[i] * B2[i] % mod; DIT(B3);
			for (i = n >> 1; i < n; ++i) b[i] = B3[n - i] * (mod - iv) % mod;
		}
	}

	// Division and Modulo Operation (use five buffers)
	void div_mod(int A, int B, pvec a, pvec b, pvec q, pvec r) {
		if (A < B) {memcpy(r, a, (A + 1) << 2), memset(r + (A + 1), 0, (B - A) << 2); return;}
		int Q = A - B, i, l_ = Q ? lg2(Q) + 1 : 0; NTT_init(l_);
		for (i = 0; i <= Q && i <= B; ++i) B4[i] = b[B - i];
		memset(B4 + i, 0, (n - i) << 2), inv(i = Q + 1, B4, B5);

		std::reverse_copy(a + B, a + (A + 1), B4), NTT_init(++l_),
		memset(B4 + i, 0, (n - i) << 2), memset(B5 + i, 0, (n - i) << 2),
		mul(2 * Q, B4, B5, q), std::reverse(q, q + (Q + 1)),
		memset(q + i, 0, (n - i) << 2);

		if (!B) return;
		NTT_init(lg2(2 * B - 1) + 1);
		for (i = 0; i <= Q && i < B; ++i) B2[i] = b[i], B3[i] = q[i];
		memset(B2 + i, 0, (n - i) << 2), memset(B3 + i, 0, (n - i) << 2),
		mul(2 * (B - 1), B2, B3, r), memset(r + i, 0, (n - i) << 2);
		for (i = 0; i < B; ++i) reduce(r[i] = a[i] - r[i]);
	}

	// Multiplication with std::vector (use two buffers, 3-dft)
	void mul(vector &a, vector &b, vector &ret) {
		int A = a.size() - 1, B = b.size() - 1;
		if (!(A || B)) {ret.EB((u64)a[0] * b[0] % mod); return;}
		NTT_init(lg2(A + B) + 1),
		memcpy(B1, a.data(), (A + 1) << 2), memset(B1 + (A + 1), 0, (n - A - 1) << 2),
		memcpy(B2, b.data(), (B + 1) << 2), memset(B2 + (B + 1), 0, (n - B - 1) << 2),
		DNTT(B1), DNTT(B2);
		for (int i = 0; i < n; ++i) B1[i] = (u64)B1[i] * B2[i] % mod;
		IDNTT(B1), ret.assign(B1, B1 + (A + B + 1));
	}

	// Differential
	void diff(int deg, pvec a, pvec b) {for (int i = 1; i <= deg; ++i) b[i - 1] = (u64)a[i] * i % mod;}

	// Integral
	void intg(int deg, pvec a, pvec b, int constant = 0) {for (int i = deg; i; --i) b[i] = (u64)a[i - 1] * ::inv[i] % mod; *b = constant;}

	// f'[x] / f[x] (use four buffers, 6.5-dft)
	void dif_quo(int deg, pvec a, pvec b) {
		assert(*a);
		if (deg <= 1) {*b = PowerMod(*a, mod - 2, a[1]); return;}

		int i, len = lg2(deg - 1);
		inv((deg + 1) / 2, a, B4), NTT_init(len + 1),
		memset(B4 + (n >> 1), 0, n << 1), DIF(B4, B2),

		diff(deg, a, B1), memcpy(B3, B1, n << 1),
		memset(B3 + (n >> 1), 0, n << 1), DIF(B3);

		for (i = 0; i < n; ++i) B3[i] = (u64)B3[i] * B2[i] % mod;
		DIT(B3, b), *b = *b * iv % mod;
		for (i = 1; i < n >> 1; ++i) b[i] = b[n - i] * iv % mod;
		memset(b + (n >> 1), 0, n << 1);

		DIF(b, B4), DIF(a, B3);
		for (i = 0; i < n; ++i) B3[i] = (u64)B3[i] * B4[i] % mod; DIT(B3);
		for (i = n >> 1; i < n; ++i) B3[i] = (B3[n - i] * iv + mod - B1[i]) % mod;

		memset(B3, 0, n << 1), DIF(B3);
		for (i = 0; i < n; ++i) B3[i] = (u64)B3[i] * B2[i] % mod; DIT(B3);
		for (i = n >> 1; i < n; ++i) b[i] = B3[n - i] * (mod - iv) % mod;
	}

	// Logarithm (use DifQuo)
	inline void ln(int deg, pvec a, pvec b) {assert(*a == 1), --deg ? (dif_quo(deg, a, b), intg(deg, b, b)) : void(*b = 0);}

	// Exponentiation (use six buffers, 12-dft)
	// WARNING : this implementation of exponentiation is SLOWER than the CDQ_NTT ver.
	void exp(int deg, pvec a, pvec b) {
		int i, len; pvec c = B6; assert(!*a);
		if (*b = 1, deg <= 1) return;
		if (b[1] = a[1], deg == 2) return;

		memset(b + 2, 0, i = 8 << lg2(deg - 1)), memset(c, 0, i), memset(B1, 0, i),
		*c = 1, neg(c[1] = b[1]);

		for (len = 1; 1 << len < deg; ++len) {
			NTT_init(len + 1);

			DIF(c, B2), DIF(b, B3);
			for (i = 0; i < n; ++i) B4[i] = (u64)B3[i] * B2[i] % mod; DIT(B4);
			for (i = n >> 1; i < n; ++i) B4[i] = B4[n - i] * iv % mod;

			memset(B4, 0, n << 1), DIF(B4);
			for (i = 0; i < n; ++i) B4[i] = (u64)B4[i] * B2[i] % mod; DIT(B4);
			for (i = n >> 1; i < n; ++i) B4[i] = B4[n - i] * (mod - iv) % mod;

			memcpy(B4, c, n << 1), DIF(B4);
			diff(n >> 1, b, B1), DIF(B1, B5);
			for (i = 0; i < n; ++i) B4[i] = (u64)B4[i] * B5[i] % mod; DIT(B4);
			for (i = n >> 1; i < n; ++i) reduce(B5[i] = (a[i] + B4[n - i + 1] * (mod - iv) % mod * ::inv[i]) % mod);

			memset(B5, 0, n << 1), DIF(B5);
			for (i = 0; i < n; ++i) B5[i] = (u64)B5[i] * B3[i] % mod; DIT(B5);
			for (i = n >> 1; i < n; ++i) b[i] = B5[n - i] * iv % mod;

			if (2 << len >= deg) return;

			DIF(b, B3);
			for (i = 0; i < n; ++i) B3[i] = (u64)B3[i] * B2[i] % mod; DIT(B3);
			for (i = n >> 1; i < n; ++i) B3[i] = B3[n - i] * iv % mod;

			memset(B3, 0, n << 1), DIF(B3);
			for (i = 0; i < n; ++i) B3[i] = (u64)B3[i] * B2[i] % mod; DIT(B3);
			for (i = n >> 1; i < n; ++i) c[i] = B3[n - i] * (mod - iv) % mod;
		}
	}
}

namespace CDQ_NTT {
	using namespace poly_base;

	int lim;
	vec f, g, C1;
	int fn[N * 2], gn[N * 2];

	inline void register_g(pvec g) {for (int i = 1; 1 << (i - 1) <= lim; ++i) NTT_init(i), DIF(g, gn + (1 << i));}

	// Standard CDQ-NTT Algorithm, type `UK`
	void solve(int L, int w) {
		int i, R = L + (1 << w), M;
		if (!w) {
			// something depend on problem
			return;
		}
		solve(L, w - 1);
		if ((M = (1 << (w - 1)) + L) > lim) return;
		NTT_init(w);
		pvec ga = gn + (1 << w);
		memcpy(C1, f + L, 2 << w), memset(C1 + (1 << (w - 1)), 0, 2 << w),
		DIF(C1);
		for (i = 0; i < n; ++i) C1[i] = (u64)C1[i] * ga[i] % mod;
		DIT(C1);
		for (i = M; i < R; ++i) f[i] = (f[i] + C1[n - (i - L)] * iv) % mod;
		solve(M, w - 1);
	}

	void solve(int L, int w) {
		int i, R = L + (1 << w), M;
		if (!w) {
			// something depend on problem
			return;
		}
		solve(L, w - 1);
		if ((M = (1 << (w - 1)) + L) > lim) return;
		NTT_init(w);
		if (L) {
			pvec fa = fn + (1 << w), ga = gn + (1 << w);
			memcpy(C1, f + L, 2 << w), memset(C1 + (1 << (w - 1)), 0, 2 << w),
			memcpy(C2, g + L, 2 << w), memset(C2 + (1 << (w - 1)), 0, 2 << w),
			DIF(C1), DIF(C2);
			for (i = 0; i < n; ++i) C1[i] = ((u64)C1[i] * ga[i] + (u64)C2[i] * fa[i]) % mod;
			DIT(C1);
			for (i = M; i < R; ++i) f[i] = (f[i] + C1[n - (i - L)] * iv) % mod;
		} else {
			memcpy(C1, f, 2 << w), memset(C1 + M, 0, 2 << w),
			memcpy(C2, g, 2 << w), memset(C2 + M, 0, 2 << w),
			DIF(C1), DIF(C2),
			memcpy(fn + (1 << (w - 1)), C1, 2 << w),
			memcpy(gn + (1 << (w - 1)), C2, 2 << w);
			for (i = 0; i < n; ++i) C1[i] = (u64)C1[i] * C2[i] % mod;
			DIT(C1);
			for (i = M; i < R; ++i) f[i] = (f[i] + C1[n - i] * iv) % mod;
		}
		solve(M, w - 1);
	}
}

namespace poly_evaluation {
	using namespace poly_base;

	int cnt = 0, lc[N], rc[N];
	vec Prd_, E1, E2, E3;
	vector g[N], tmp_;

	int solve(int L, int R) {
		if (L + 1 == R) return L;
		int M = (L + R) / 2, id = cnt++, lp = solve(L, M), rp = solve(M, R);
		return poly::mul(g[lp], g[rp], g[id]), lc[id] = lp, rc[id] = rp, id;
	}

	void recursion(int id, int L, int R, const vector &poly) {
		if (L + 1 == R) return tmp_.EB(poly.back());
		int i, n = poly.size() - 1, M = (L + R) / 2, lp = lc[id], rp = rc[id],
			dl = min(n, g[lp].size() - 1), dr = min(n, g[rp].size() - 1);
		if (L + 2 == R) return
			tmp_.EB((poly[n] + (u64)poly[n - 1] * g[rp].back()) % mod),
			tmp_.EB((poly[n] + (u64)poly[n - 1] * g[lp].back()) % mod);

		vector ly, ry; ly.reserve(dl + 1), ry.reserve(dr + 1);
		NTT_init(lg2(dl + dr) + 1);
		memcpy(E1, poly.data(), (n + 1) << 2), DIF(E1, E2), memset(E1, 0, (n + 1) << 2);

		memcpy(E1, g[rp].data(), (dr + 1) << 2), DIF(E1, E3), memset(E1, 0, (dr + 1) << 2);
		for (i = 0; i < poly::n; ++i) E3[i] = (u64)E3[i] * E2[i] % mod;
		DIT(E3), std::reverse(E3 + 1, E3 + poly::n);
		for (i = n - dl; i <= n; ++i) ly.EB(E3[i] * iv % mod);

		memcpy(E1, g[lp].data(), (dl + 1) << 2), DIF(E1, E3), memset(E1, 0, (dl + 1) << 2);
		for (i = 0; i < poly::n; ++i) E3[i] = (u64)E3[i] * E2[i] % mod;
		DIT(E3), std::reverse(E3 + 1, E3 + poly::n);
		for (i = n - dr; i <= n; ++i) ry.EB(E3[i] * iv % mod);

		recursion(lp, L, M, ly), recursion(rp, M, R, ry);
	}

	vector emain(int n, pvec f, const vector &pts) {
		int i, id, m = pts.size(), q;
		if (!m) return vector();
		if (!n) return vector(m, *f);
		for (i = 0; i < m; ++i) g[i].clear(), g[i].EB(1), g[i].EB(neg(q = pts[i]));

		id = solve(0, cnt = m), memcpy(Prd_, g[id].data(), (m + 1) << 2);
		poly::inv(n + 1, Prd_, E2), memset(Prd_, 0, (m + 1) << 2);
		if (n > 0) memset(E2 + (n + 1), 0, (poly::n - n - 1) << 2);

		std::reverse_copy(f, f + (n + 1), E1), poly::mul(2 * n, E1, E2, E3),
		memset(E1, 0, (n + 1) << 2), memset(E2, 0, (n + 1) << 2);

		return tmp_.clear(), tmp_.reserve(m), recursion(id, 0, m, vector(E3 + max(n - m + 1, 0), E3 + (n + 1))), tmp_;
	}
}

namespace poly_interpolation {
	using namespace poly_evaluation;

	vec I1, I2, I3, I4, I5;
	pvec iresult;

	void irecursion(int id, int L, int R) {
		if (L + 1 == R) return;
		int i, M = (L + R) / 2, lp = lc[id], rp = rc[id];
		irecursion(lp, L, M), irecursion(rp, M, R);
		NTT_init(lg2(R - L - 1) + 1);
		memcpy(I1, iresult + L, (M - L) << 2), memset(I1 + (M - L), 0, (poly::n - (M - L)) << 2);
		memcpy(I2, iresult + M, (R - M) << 2), memset(I2 + (R - M), 0, (poly::n - (R - M)) << 2);
		memcpy(I3, g[lp].data(), (M - L + 1) << 2), memset(I3 + (M - L + 1), 0, (poly::n - (M - L + 1)) << 2);
		memcpy(I4, g[rp].data(), (R - M + 1) << 2), memset(I4 + (R - M + 1), 0, (poly::n - (R - M + 1)) << 2);
		DIF(I1), DIF(I2), DIF(I3), DIF(I4);
		for (i = 0; i < n; ++i) I1[i] = ((u64)I1[i] * I4[i] + (u64)I2[i] * I3[i]) % mod;
		DIT(I1), std::reverse(I1 + 1, I1 + n);
		for (i = 0; i < R - L; ++i) iresult[L + i] = I1[i] * iv % mod;
	}

	void imain(int n, pr *pts, pvec ret) {
		int i, id, q;
		assert(n > 0);
		for (i = 0; i < n; ++i) g[i].clear(), g[i].EB(1), g[i].EB(neg(q = pts[i].first));

		id = solve(0, cnt = n), memcpy(Prd_, g[id].data(), (n + 1) << 2);
		poly::inv(n, Prd_, E2), memset(Prd_, 0, n << 2);
		if (n > 1) memset(E2 + n, 0, (poly::n - n) << 2);

		for (i = 0; i < n; ++i) E1[i] = u64(n - i) * g[id][i] % mod;
		poly::mul(2 * (n - 1), E1, E2, E3),
		memset(E1, 0, n << 2), memset(E2, 0, n << 2);

		tmp_.clear(), tmp_.reserve(n), recursion(id, 0, n, vector(E3, E3 + n));
		for (i = 0; i < n; ++i) ret[i] = PowerMod(tmp_[i], mod - 2, pts[i].second);
		iresult = ret, irecursion(id, 0, n), std::reverse(ret, ret + n);
	}
}

namespace linear_recur {
	using namespace poly_base;

	int ld;
	vec f, g, L1, L2, L3;

	void __builtin_divmod(int d) {
		int i; NTT_init(ld);
		std::reverse_copy(f + d, f + 2 * d, L3);
		memset(L3 + d, 0, (n - d) << 2);
		DNTT(L3);
		for (i = 0; i < n; ++i) L3[i] = (u64)L3[i] * L2[i] % mod;
		IDNTT(L3);
		std::reverse(L3, L3 + d), memset(L3 + d, 0, (n - d) << 2);
		DNTT(L3);
		for (i = 0; i < n; ++i) L3[i] = (u64)L3[i] * L1[i] % mod;
		IDNTT(L3);
		for (i = 0; i < d; ++i) reduce(f[i] -= L3[i]);
		memset(f + d, 0, d << 2);
	}

	void solve(int Q, int d) {
		int i, df = 1, z = lg2(Q); bool alive = false;
		assert(d > 0), ld = lg2(2 * d - 1) + 1;
		if (d == 1) alive = true, *f = 1;

		NTT_init(ld),
		std::reverse_copy(g + 1, g + (d + 1), L1),
		memset(L1 + d, 0, (n - d) << 2),
		poly::inv(d, L1, L2);

		NTT_init(ld),
		memset(L2 + d, 0, (n - d) << 2),
		memcpy(L1, g, d << 2),
		memset(L1 + d, 0, (n - d) << 2),
		DIF(L1), DIF(L2);

		for (; --z >= 0; ) {
			df = df << 1 | (Q >> z & 1);
			if (df < d) continue;
			if (!alive) alive = true, f[df >> 1] = 1;
			NTT_init(ld), DNTT(f);
			for (i = 0; i < n; ++i) f[i] = (u64)f[i] * f[i] % mod;
			IDNTT(f);
			if (df & 1) std::copy_backward(f, f + (2 * d - 1), f + 2 * d), *f = 0;
			__builtin_divmod(d);
		}
	}
}

namespace linear_recur_single {
	using namespace poly;

	int ld;
	vec f, g, L1, L2, L3, L4;

	int get(int Q, int d) {
		if (Q < d) return f[Q];

		int i, j, ret = 0, halves = 0;
		assert(d > 0), ld = lg2(2 * d - 1) + 1,

		mul(2 * d - 1, f, g, L1), memset(L1 + d, 0, (n - d) << 2),
		memcpy(L2, g, (d + 1) << 2), memset(L2 + (d + 1), 0, (n - d - 1) << 2);

		for (; Q >= d; Q >>= 1) {
			NTT_init(ld), DIF(L1, L4), DIF(L2, L3),
			NTT_init(ld - 1), ++halves;
			if (Q & 1) {
				for (j = i = 0; i < n; ++i, j += 2)
					L4[i] = ((u64)L4[j] * L3[j + 1] + u64(mod - L3[j]) * L4[j + 1]) % mod;
				for (i = 0; i < n; ++i) L4[i] = (u64)L4[i] * w2[i] % mod;
				IDNTT(L4), L4[n] = *L4, memcpy(L1, L4 + 1, d << 2);
			} else {
				for (j = i = 0; i < n; ++i, j += 2)
					L4[i] = ((u64)L4[j] * L3[j + 1] + (u64)L3[j] * L4[j + 1]) % mod;
				IDNTT(L4, L1);
			}
			for (j = i = 0; i < n; ++i, j += 2) L3[i] = (u64)L3[j] * L3[j + 1] % mod;
			IDNTT(L3, L2);
			if (*L2 != 1) L2[d] = (*L2 ? *L2 : mod) - 1, *L2 = 1;
		}

		poly::inv(Q + 1, L2, L3);
		for (i = 0; i <= Q; ++i) ret = (ret + (u64)L1[i] * L3[Q - i]) % mod;
		return PowerMod(iv2, halves, ret);
	}
}

namespace miscellaneous {
	using namespace poly_base;

	vec M1, M2, M3, M4;
	vec Ma[200];

	// calculate b(x) = a(x + c), poly_base::init need 'false'
	void translate(int deg, pvec a, int c, pvec b) {
		if (!deg) {*b = *a; return;}
		int i; u64 pc = 1; NTT_init(lg2(deg) + 2);
		for (i = 0; i <= deg; ++i) M1[i] = (u64)a[i] * fact[i] % mod;
		memset(M1 + i, 0, (n - i) << 2);
		for (i = deg; i >= 0; --i, pc = pc * c % mod) M2[i] = finv[deg - i] * pc % mod;
		DIF(M1), DIF(M2);
		for (i = 0; i < n; ++i) M1[i] = (u64)M1[i] * M2[i] % mod;
		DIT(M1);
		for (i = 0; i <= deg; ++i) b[i] = M1[n - deg - i] * iv % mod * finv[i] % mod;
	}

	// compute a(1), a(z), a(z^2), ..., a(z^(len-1)) to b
	void czt(int deg, int len, int z, pvec a, pvec b) {
		int i; u64 iz = PowerMod(z, mod - 2);
		if (!len) return;
		if (!deg || !z) return std::fill(b, b + len, *a);
		NTT_init(lg2(deg + len - 1) + 1);
		M1[1] = *M1 = M3[1] = *M3 = 1;
		for (i = 2; i < deg + len; ++i) M1[i] = M1[i - 1] * (u64)z % mod, M3[i] = M3[i - 1] * iz % mod;
		for (i = 2; i < deg + len; ++i) M1[i] = (u64)M1[i] * M1[i - 1] % mod, M3[i] = (u64)M3[i] * M3[i - 1] % mod;
		memset(M1 + i, 0, (n - i) << 2);
		for (i = 0; i <= deg; ++i) M2[deg - i] = (u64)M3[i] * a[i] % mod;
		memset(M2 + i, 0, (n - i) << 2),
		DIF(M1), DIF(M2);
		for (i = 0; i < n; ++i) M1[i] = (u64)M1[i] * M2[i] % mod;
		DIT(M1);
		for (i = 0; i < len; ++i) b[i] = M1[n - deg - i] * iv % mod * M3[i] % mod;
	}

	inline int __builtin_inner(int len, pvec a, pvec b) {
		int i = 0; u64 ans = 0;
		#define term(k) (u64)a[i + k] * b[i + k]
		for (; i + 15 < len; i += 16) {
			ans = (ans + term(0) + term(1) + term(2) + term(3)
					   + term(4) + term(5) + term(6) + term(7)
					   + term(8) + term(9) + term(10) + term(11)
					   + term(12) + term(13) + term(14) + term(15)
				) % mod;
		}
		#undef term
		for (; i < len; ++i) ans += (u64)a[i] * b[i];
		return ans % mod;
	}

	// compositional inverse
	void cinv(int deg, pvec a, pvec b) {
		assert(!*a), *b = 0;
		if (--deg <= 0) return;
		int i, j, k, d, B, *c;

		NTT_init(lg2(2 * deg - 1)),
		memcpy(M1, a + 1, deg << 2), memset(M1 + deg, 0, (n - deg) << 2),
		poly::inv(deg, M1, M2),

		NTT_init(lg2(2 * deg - 1) + 1),
		memset(M2 + deg, 0, (n - deg) << 2),

		DIF(M2, M1), B = sqrt(deg),
		memset(*Ma, 0, n << 2), **Ma = 1;
		for (j = 1; j <= B; ++j) {
			DIF(Ma[j - 1], M3), c = Ma[j];
			for (i = 0; i < n; ++i) M3[i] = (u64)M3[i] * M1[i] % mod;
			DIT(M3), *c = *M3 * iv % mod;
			for (i = 1; i < deg; ++i) c[i] = M3[n - i] * iv % mod;
			memset(c + i, 0, (n - i) << 2);
		}
		DIF(Ma[B], M1), memset(M3, 0, n << 2), *M3 = M4[deg - 1] = 1;

		for (d = 0; ; ) {
			for (k = 0; k < B && d < deg; ++k, ++d)
				b[d + 1] = (u64)__builtin_inner(d + 1, Ma[k + 1], M4 + (deg - d - 1)) * inv[d + 1] % mod;
			if (d >= deg) break;
			DIF(M3);
			for (i = 0; i < n; ++i) M3[i] = (u64)M3[i] * M1[i] % mod;
			DIT(M3), *M3 = *M3 * iv % mod;
			for (i = 1; i < deg; ++i) M3[i] = M3[n - i] * iv % mod;
			memset(M3 + i, 0, (n - i) << 2),
			std::reverse_copy(M3, M3 + deg, M4);
		}
	}
}


  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91
  • 92
  • 93
  • 94
  • 95
  • 96
  • 97
  • 98
  • 99
  • 100
  • 101
  • 102
  • 103
  • 104
  • 105
  • 106
  • 107
  • 108
  • 109
  • 110
  • 111
  • 112
  • 113
  • 114
  • 115
  • 116
  • 117
  • 118
  • 119
  • 120
  • 121
  • 122
  • 123
  • 124
  • 125
  • 126
  • 127
  • 128
  • 129
  • 130
  • 131
  • 132
  • 133
  • 134
  • 135
  • 136
  • 137
  • 138
  • 139
  • 140
  • 141
  • 142
  • 143
  • 144
  • 145
  • 146
  • 147
  • 148
  • 149
  • 150
  • 151
  • 152
  • 153
  • 154
  • 155
  • 156
  • 157
  • 158
  • 159
  • 160
  • 161
  • 162
  • 163
  • 164
  • 165
  • 166
  • 167
  • 168
  • 169
  • 170
  • 171
  • 172
  • 173
  • 174
  • 175
  • 176
  • 177
  • 178
  • 179
  • 180
  • 181
  • 182
  • 183
  • 184
  • 185
  • 186
  • 187
  • 188
  • 189
  • 190
  • 191
  • 192
  • 193
  • 194
  • 195
  • 196
  • 197
  • 198
  • 199
  • 200
  • 201
  • 202
  • 203
  • 204
  • 205
  • 206
  • 207
  • 208
  • 209
  • 210
  • 211
  • 212
  • 213
  • 214
  • 215
  • 216
  • 217
  • 218
  • 219
  • 220
  • 221
  • 222
  • 223
  • 224
  • 225
  • 226
  • 227
  • 228
  • 229
  • 230
  • 231
  • 232
  • 233
  • 234
  • 235
  • 236
  • 237
  • 238
  • 239
  • 240
  • 241
  • 242
  • 243
  • 244
  • 245
  • 246
  • 247
  • 248
  • 249
  • 250
  • 251
  • 252
  • 253
  • 254
  • 255
  • 256
  • 257
  • 258
  • 259
  • 260
  • 261
  • 262
  • 263
  • 264
  • 265
  • 266
  • 267
  • 268
  • 269
  • 270
  • 271
  • 272
  • 273
  • 274
  • 275
  • 276
  • 277
  • 278
  • 279
  • 280
  • 281
  • 282
  • 283
  • 284
  • 285
  • 286
  • 287
  • 288
  • 289
  • 290
  • 291
  • 292
  • 293
  • 294
  • 295
  • 296
  • 297
  • 298
  • 299
  • 300
  • 301
  • 302
  • 303
  • 304
  • 305
  • 306
  • 307
  • 308
  • 309
  • 310
  • 311
  • 312
  • 313
  • 314
  • 315
  • 316
  • 317
  • 318
  • 319
  • 320
  • 321
  • 322
  • 323
  • 324
  • 325
  • 326
  • 327
  • 328
  • 329
  • 330
  • 331
  • 332
  • 333
  • 334
  • 335
  • 336
  • 337
  • 338
  • 339
  • 340
  • 341
  • 342
  • 343
  • 344
  • 345
  • 346
  • 347
  • 348
  • 349
  • 350
  • 351
  • 352
  • 353
  • 354
  • 355
  • 356
  • 357
  • 358
  • 359
  • 360
  • 361
  • 362
  • 363
  • 364
  • 365
  • 366
  • 367
  • 368
  • 369
  • 370
  • 371
  • 372
  • 373
  • 374
  • 375
  • 376
  • 377
  • 378
  • 379
  • 380
  • 381
  • 382
  • 383
  • 384
  • 385
  • 386
  • 387
  • 388
  • 389
  • 390
  • 391
  • 392
  • 393
  • 394
  • 395
  • 396
  • 397
  • 398
  • 399
  • 400
  • 401
  • 402
  • 403
  • 404
  • 405
  • 406
  • 407
  • 408
  • 409
  • 410
  • 411
  • 412
  • 413
  • 414
  • 415
  • 416
  • 417
  • 418
  • 419
  • 420
  • 421
  • 422
  • 423
  • 424
  • 425
  • 426
  • 427
  • 428
  • 429
  • 430
  • 431
  • 432
  • 433
  • 434
  • 435
  • 436
  • 437
  • 438
  • 439
  • 440
  • 441
  • 442
  • 443
  • 444
  • 445
  • 446
  • 447
  • 448
  • 449
  • 450
  • 451
  • 452
  • 453
  • 454
  • 455
  • 456
  • 457
  • 458
  • 459
  • 460
  • 461
  • 462
  • 463
  • 464
  • 465
  • 466
  • 467
  • 468
  • 469
  • 470
  • 471
  • 472
  • 473
  • 474
  • 475
  • 476
  • 477
  • 478
  • 479
  • 480
  • 481
  • 482
  • 483
  • 484
  • 485
  • 486
  • 487
  • 488
  • 489
  • 490
  • 491
  • 492
  • 493
  • 494
  • 495
  • 496
  • 497
  • 498
  • 499
  • 500
  • 501
  • 502
  • 503
  • 504
  • 505
  • 506
  • 507
  • 508
  • 509
  • 510
  • 511
  • 512
  • 513
  • 514
  • 515
  • 516
  • 517
  • 518
  • 519
  • 520
  • 521
  • 522
  • 523
  • 524
  • 525
  • 526
  • 527
  • 528
  • 529
  • 530
  • 531
  • 532
  • 533
  • 534
  • 535
  • 536
  • 537
  • 538
  • 539
  • 540
  • 541
  • 542
  • 543
  • 544
  • 545
  • 546
  • 547

参考资料

  1. 学习笔记：多项式全家桶 By NaCly_Fish %鱼
  2. 多项式全家桶
  3. FFT什么的 By YWW
  4. 【Learning】 多项式全家桶（加减乘除逆元ln exp）ByRogerDTZ
  5. OI Wiki _多项式 OI Wiki yyds！
  6. [学习笔记]多项式相关运算2(求导/积分,牛顿迭代,开根,求ln/exp,快速幂,拉格朗日反演)
  7. 浅谈多项式 by CYJian
声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/爱喝兽奶帝天荒/article/detail/948299
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号