FFT
坑
NTT
将(FFT)中的单位复数根改成原根即可。
卡常版NTT模版
struct Mul
{ int Len;
int wn[N], Lim;
int rev[N];
inline void getReverse(int * a)
{
static int rev[N];
rev[0] = 0;
for (register int i = 0; i < Len; i++)
{
rev[i] = (rev[i>>1] >> 1) | ((i&1) ? (Len >> 1) : 0);
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
}
void NTT(int * a, int opt)
{
getReverse(a);
for (register int i = 0; i < Len; i += 2)
{
int x = a[i], y = 1LL * wn[0] * a[i + 1] % mod;
a[i] = Plus(x, y), a[i + 1] = Minus(x, y);
}
for (register int i = 0; i < Len; i += 4)
{
int x1 = a[i], y1 = 1LL * wn[0] * a[i + 2] % mod;
int x2 = a[i + 1], y2 = 1LL * wn[Lim / 2] * a[i + 3] % mod;
a[i] = Plus(x1, y1), a[i + 2] = Minus(x1, y1);
a[i + 1] = Plus(x2, y2), a[i + 3] = Minus(x2, y2);
}
for (register int i = 8; i <= Len; i <<= 1)
{
int M = i >> 1;
for (register int j = 0; j < Len; j += i)
{
for (register int k = 0; k < M; k += 4)
{
int x1 = a[j + k], y1 = 1LL * wn[Lim / M * k] * a[j + k + M] % mod;
int x2 = a[j + k + 1], y2 = 1LL * wn[Lim / M * (k+1)] * a[j + k + M + 1] % mod;
int x3 = a[j + k + 2], y3 = 1LL * wn[Lim / M * (k+2)] * a[j + k + M + 2] % mod;
int x4 = a[j + k + 3], y4 = 1LL * wn[Lim / M * (k+3)] * a[j + k + M + 3] % mod;
a[j + k] = Plus(x1, y1), a[j + k + M] = Minus(x1, y1);
a[j + k + 1] = Plus(x2, y2), a[j + k + M + 1] = Minus(x2, y2);
a[j + k + 2] = Plus(x3, y3), a[j + k + M + 2] = Minus(x3, y3);
a[j + k + 3] = Plus(x4, y4), a[j + k + M + 3] = Minus(x4, y4);
}
}
}
if (opt == -1)
{
for (register int i = 0; i < Len; i++)
a[i] = 1LL * a[i] * inv[Len] % mod;
for(register int i = 1; i < Len; i++)
if(i < Len-i) swap(a[i], a[Len-i]);
}
}
inline void init()
{
Lim = Len;
wn[0] = 1; wn[1] = power(7, (mod-1) / (Len << 1), mod);
for (int i = 2; i < (Len << 1); i++) wn[i] = 1LL * wn[i-1] * wn[1] % mod;
}
inline int getLen(int l)
{
for (Len = 1; Len <= l; Len <<= 1);
return Len;
}
};
多项式求逆
也就是求
[A(x)B(x) = 1 pmod{x^n}
]
假设我们已知
[A(x)G(x) = 1 pmod{x^{lceil{n over 2}
ceil}}
]
两式相减得
[A(x)(B(x) - G(x)) = 0 pmod{x^{lceil{n over 2}
ceil}}
]
[B(x) - G(x) = 0 pmod{x^{lceil{n over 2}
ceil}}
]
[B^2(x) + G^2(x) - 2B(x)G(x) = 0 pmod{x^n}
]
两边同时乘上(A(x))
[A(x)B^2(x) + A(x)G^2(x) = 2A(x)B(x)G(x)
]
因为(A(x)B(x) = 1), 所以可以得到
[B(x) + A(x)G^2(x) = 2G(x)
]
[B(x) = 2G(x) - A(x)G^2(x)
]
就这样解决了
时间复杂度为(T(n) = T(frac{n}{2}) +O(n log n))
代码
inline void getInv(int * A, int * B, int len)
{
static int tmp1[N], tmp2[N];
B[0] = power(A[0], mod-2, mod);
for (register int k = 2; k <= len; k <<= 1)
{
int Len = k << 1;
cpy(tmp1, A, k, Len);
cpy(tmp2, B, (k >> 1), Len);
Calc.Len = Len;
Calc.NTT(tmp1, 1);
Calc.NTT(tmp2, 1);
for (int i = 0; i < Len; i++)
tmp1[i] = Minus(Plus(tmp2[i], tmp2[i]), 1LL * tmp1[i] * tmp2[i] % mod * tmp2[i] % mod);
Calc.NTT(tmp1, -1);
cpy(B, tmp1, k, Len);
}
return;
}
多项式开根
也就是求
[B^2(x) = A(x) pmod{x^n}
]
假设我们已知
[G^2(x) = A(x) pmod{x^{lceil{n over 2}
ceil}}
]
两式相减得
[B^2(x) - G^2(x) = 0 pmod{x^{lceil{n over 2}
ceil}}
]
然后平方一下
[B^4(x) + G^4(x) - 2B^2(x)G^2(x) = 0 pmod{x^n}
]
[B^4(x) + G^4(x) = 2B^2(x)G^2(x) pmod{x^n}
]
配一下方
[B^4(x) + G^4(x) + 2B^2(x)G^2(x) = 4B^2(x)G^2(x) pmod{x^n}
]
[(B^2(x) + G^2(x))^2 = (2B(x)G(x))^2 pmod{x^n}
]
[B^2(x) + G^2(x) = 2B(x)G(x) pmod{x^n}
]
因为(B^2(x) = A(x)), 所以可以得到
[B(x) = {A(x) + G^2(x) over {2G(x)}}
]
为了方便实现，我们常把它化成
[B(x) = {A(x) over 2G(x)} + {G(x) over 2}
]
然后就解决了
时间复杂度为(T(n) = T(frac{n}{2}) + O(n log n) = O(n log n))
多项式求导
模拟即可
代码
inline void getDeri(int * a, int len)
{ for (register int i = 0; i < len; i++)
a[i] = 1LL * a[i+1] * (i+1) % mod;
}
多项式积分
模拟即可
代码
inline void getInte(int * a, int len)
{ for (int i = len-1; i >= 1; i--)
a[i] = 1LL * a[i-1] * inv[i] % mod;
a[0] = 0;
}
多项式求(ln)
已知多项式(A(x)), 求(B(x) = ln (A(x)))
根据链式法则， 我们可以得到
[B'(x) = frac{mathbb{d}(ln(A(x)))}{mathbb{d}A(x)} frac{mathbb{d}A(x)}{mathbb{d}x} = frac{A'(x)}{A(x)}
]
所以
[B(x) = int frac{A'(x)}{A(x)} mathbb{d}x
]
代码
void getLn(int * A, int len)
{
static int tmp1[N], tmp2[N], tmp3[N];
int Len = len << 1;
cpy(tmp1, A, len, Len);
cpy(tmp2, A, len, Len);
getDeri(tmp1, len);
getInv(tmp2, tmp3, len);
Calc.Len = Len;
Calc.NTT(tmp1, 1);
Calc.NTT(tmp3, 1);
for (int i = 0; i < Len; i++)
tmp1[i] = 1LL * tmp1[i] * tmp3[i] % mod;
Calc.NTT(tmp1, -1);
memset(tmp1 + len, 0, 4 * (Len - len));
getInte(tmp1, len);
cpy(A, tmp1, len, Len);
}
时间复杂度还是(T(n) = T(frac{n}{2}) + O(n log n) = O(n log n))
牛顿迭代
我们要求(f(x) = 0)的根
那么可以使用泰勒公式来近似(f(x) = 0)的根
我们设当前的近似值是(x_{n})， 我们想要得到的近似值是(x_{n+1})
截取泰勒公式的线性部分： $$f(x_{n}) + f'(x_{n})(x_{n+1}-x_n) = 0$$
解方程得： $$x_{n+1} = x_n - frac{f(x_{n})}{f'(x_{n})}$$
这样不断迭代
我们可以用这种方法来解多项式方程。
多项式(exp)
我们要求的是 (exp(A(x))), 这相当于解一个方程： (ln(exp(A(x))) - A(x) = 0)
可以直接套用牛顿迭代法求解。