多项式全家桶（但还没学明白，这篇笔记也还没完工）
Crescent_Rose_ · · 算法·理论
XK 的课件（其中可能有一些小错误）
stO XK Orz
递归 FFT
#include <bits/stdc++.h>
//#define int long long
#define Comp complex < double > //
using namespace std;
namespace IO
{
int read()
{
int f = 1, x = 0;
char c = getchar();
while(c < '0' || c > '9'){
if(c == '-') f = - 1;
c = getchar();
}
while(c >= '0' && c <= '9'){
x = x * 10 + c - '0';
c = getchar();
}
return f * x;
}
void write(int x)
{
if(x < 0){
putchar('-');
x = - x;
}
if(x > 9) write(x / 10);
putchar(x % 10 + '0');
}
}
using namespace IO;
// Fast IO sometimes -> Slow IO
const int N = (1 << 21);
const double PI = acos(- 1); // 不知道为什么
Comp ff[N + 1], g[N + 1];
void DFT(Comp * f, int n, int rev)
{
if(n == 1) return;
Comp f1[n / 2], f2[n / 2];
for(int i = 0; i < n / 2; i ++){
f1[i] = f[i * 2];
f2[i] = f[i * 2 + 1];
}
DFT(f1, n / 2, rev);
DFT(f2, n / 2, rev);
Comp wk(1, 0), w1(cos((double)PI * 2 / n), sin((double)PI * 2 / n) * (double)rev); //
// 为什么用下面注释掉的，不用后面没注释掉的，就不对？
// int i;
// for(i = 0; i < n / 2; i ++){
// f[i] = f1[i] + wk * f2[i];
// wk *= w1;
// }
// for(; i < n; i ++){
// f[i] = f[i - n / 2] + wk * f2[i - n / 2];
// wk *= w1;
// }
for(int i = 0; i < n / 2; i ++){
f[i] = f1[i] + wk * f2[i];
f[i + n / 2] = f1[i] - wk * f2[i]; //
wk *= w1;
}
}
void solve()
{
int n, m;
n = read(); m = read();
for(int i = 0; i <= n; i ++){
double tmp = read();
ff[i].real(tmp);
}
for(int i = 0; i <= m; i ++){
double tmp = read();
g[i].real(tmp);
}
int len = 1;
while(len < n + m + 1) len <<= 1; // (n + m) 次，一共 (n + m + 1) 项（？）
DFT(ff, len, 1);
DFT(g, len, 1);
for(int i = 0; i < len; i ++) ff[i] *= g[i];
DFT(ff, len, - 1);
for(int i = 0; i <= n + m; i ++) printf("%.0lf ", 1e-7 + ff[i].real() / (double)len); // 为什么要加一个很小的数？为了补上丢失的精度？
putchar('\n');
}
signed main()
{
int t = 1;
while(t --) solve();
return 0;
}
// 递归 FFT非递归 FFT
#include <bits/stdc++.h>
#define Complex complex < double > //
#define PI acos(-1) // ???
using namespace std;
const int L = (1 << 21); // ?
Complex f[L], g[L];
int rev[L]; // ?
void Change(Complex y[], int len)
{
for(int i = 0; i < len; ++ i){
rev[i] = rev[i >> 1] >> 1;
if(i & 1) rev[i] |= (len >> 1);
}
for(int i = 0; i < len; ++ i) if(i < rev[i]) swap(y[i], y[rev[i]]);
}
void DFT(Complex y[], int len, int Rev)
{
Change(y, len);
for(int h = 2; h <= len; h <<= 1){
Complex wn(cos(2 * PI / h), Rev * sin(2 * PI / h)); //
for(int j = 0; j < len; j += h){
Complex w(1, 0);
for(int k = j; k < j + h / 2; k ++){
Complex u = y[k], t = w * y[k + h / 2];
y[k] = u + t; y[k + h / 2] = u - t;
w *= wn;
}
}
}
}
int main()
{
int n, m; scanf("%d%d", & n, & m);
for(int i = 0; i <= n; i ++){
int tmp; scanf("%d", & tmp);
Complex t(tmp, 0);
f[i] = t;
}
for(int i = 0; i <= m; i ++){
int tmp; scanf("%d", & tmp);
Complex t(tmp, 0);
g[i] = t;
}
int len = (1 << ((int)ceil(log2((n + m + 1))))); // ?
DFT(f, len, 1);
DFT(g, len, 1);
for(int i = 0; i < len; i ++) f[i] *= g[i]; // ???
DFT(f, len, - 1);
for(int i = 0; i < len; i ++) f[i] /= len; // ???
for(int i = 0; i < n + m + 1; i ++) printf("%.0lf ", (f[i].real() + 1e-7));
return 0;
}非递归 NTT
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int L = (1 << 21), P = 998244353, G = 3; // ?
int f[L], g[L];
int rev[L]; // ?
int Power(int x, int y, int p)
{
int res = 1; // 是 1 而非 0
while(y){
if(y & 1) res = res * x % p;
x = x * x % p;
y >>= 1;
}
return res;
}
void Change(int y[], int len)
{
for(int i = 0; i < len; ++ i){
rev[i] = rev[i >> 1] >> 1;
if(i & 1) rev[i] |= (len >> 1);
}
for(int i = 0; i < len; ++ i) if(i < rev[i]) swap(y[i], y[rev[i]]);
}
void NTT(int y[], int len, int Rev)
{
Change(y, len);
for(int h = 2; h <= len; h <<= 1){
int wn = Power(G, (P - 1) / h, P); // 是 / h 而不是 / len
if(Rev == - 1) wn = Power(wn, P - 2, P); // ?
for(int j = 0; j < len; j += h){
int w = 1;
for(int k = j; k < j + h / 2; k ++){
int u = y[k], t = w * y[k + h / 2] % P;
y[k] = (u + t) % P; y[k + h / 2] = (u - t + P) % P;
w = w * wn % P;
}
}
}
}
signed main()
{
int n, m; scanf("%lld%lld", & n, & m);
for(int i = 0; i <= n; i ++) scanf("%lld", & f[i]); //
for(int i = 0; i <= m; i ++) scanf("%lld", & g[i]); //
int len = (1 << ((int)ceil(log2((n + m + 1))))); // ?
NTT(f, len, 1);
NTT(g, len, 1);
for(int i = 0; i < len; i ++) f[i] = f[i] * g[i] % P; // ???
NTT(f, len, - 1);
int inv = Power(len, P - 2, P); // 把下一句中的 inv 拿到这里算会快不少
for(int i = 0; i < n + m + 1; i ++) f[i] = f[i] * inv % P; // 要取模
for(int i = 0; i < n + m + 1; i ++) printf("%lld ", f[i]);
return 0;
}加法
略。
乘法
略。
求逆
求
sol:
牛顿迭代的一个重要式子：
每次
还要考虑一开始的
据 lr：
不确定是不是
现在来求逆。
把