$FFT前置知识及FFT\FFT分治$

FFT前置知识及FFT\FFT分治

2020年02月26日算法竞赛

ICPC 数论

首先依次介绍FFT(快速傅里叶变换)前置知识: 复数及单位根
复变函数欧拉定理单位根三大引理多项式的系数与点值表示法多项式的求值与插值(顺介绍拉格朗日插值法) 向量卷积

复数及单位根

我们知道复数可以写成,其中为实部,为虚部.而且一个复数可以在一个复平面上表示出来,而我们可以将其写成极坐标的形式而根据欧拉定理:

则有:,这里给出乘法的例子:

下面介绍单位根: 如果存在复数,则其根被称作单位根,显然单位根的模为1,其中我们定义单位根:

其中定义主n次单位根(当k==1时):

而在程序中我们一般写成三角函数形式:

这里给出复数运算模板代码:

struct Complex
{
    double r, i;
    Complex() {}//无参函数不要写东西!易超时
    Complex(double r, double i)
    :r(r), i(i) {}
}f[N], g[N];
Complex operator + (const Complex & a, const Complex & b)
{return Complex(a.r + b.r, a.i + b.i);}
Complex operator - (const Complex & a, const Complex & b)
{return Complex(a.r - b.r, a.i - b.i);}
Complex operator * (const Complex & a, const Complex & b)

下面给出单位根的三大引理

多项式的两种表示方法: 第一种非常常见的系数表示法,及,,另外我们给出另一种表示方法,点值表示法,及用n+1种不同的取值代入,组成的集合:

我们分析两种表示方法: 系数表达式: 求值:我们可以使用秦九昭算法得出, 求和:两个表达式的各项系数相加,是的求积:一个表达式相乘的各项都需要和另一个的各项相乘,是的 点值表达式: 求值:我们使用拉格朗日插值法可以在的出系数表达式并求值(插值的定义一会给出) 这里给出拉格朗日插值法的公式(n项表达式):

并有模板: P4781 【模板】拉格朗日插值下面是模板代码:

const int N = 2e3+5;
const ll mod = 998244353;
ll a[N], b[N];
ll inv(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if (b&1) res = res * a % mod;
        b >>= 1;
        a = a * a % mod;
    }
    return res % mod;
}
int main()
{
    int n; ll k;
    cin >> n >> k;
    for (int i = 1; i <= n; i++)
        scanf("%d%d", &a[i], &b[i]);
    ll sum = 0;
    for (int i = 1; i <= n; i++)
    {
        ll ans = b[i];
        ll fa = 1, fb = 1;
        for (int j = 1; j <= n; j++)
        {
            if (j == i) continue;
            fa = (fa * (k - a[j] + mod)%mod)%mod;
            fb = (fb * (a[i] - a[j] + mod)%mod)%mod;
        }
        ans = (ans * fa % mod * inv(fb, mod-2)%mod)%mod;
        sum = (sum + ans)%mod;
    }
    printf("%lld\n", sum);
}

求和:两个表达式的各项点值的

值相加,是

的求积:两个表达式的各项点值的

值相乘,是

的我们不难发现,如果我们相求两个多项式的积

(系数向量分别是

),在点值表达式的情况下将十分高效.这时新表达式

(系数向量是c)的第i项系数为:

我们则称为的卷积,记作,两个多项式相乘,就是对应的两个系数向量的卷积, 我们将一个点值表达式转化为系数表达式称为插值. 由此可见:如何提高多项式求值与插值的效率.决定着两个表达式相乘的效率!

话讲完,正式下面介绍FFT

首先看DFT(离散傅里叶变换) 我们将单位根带入到系数表达式(方便起见下式中用替代)中:

记作每次求值为所以离散傅里叶变换最终时间复杂度为这样显然是不行的,毕竟直接相乘的时间复杂度都是然后就有了FFT: 我们对于一个n项系数表达式(n为2的幂数):

可以分成:

这时我们发现一个及其有趣的事实:

而根据(1),(2)两个引理:我们最终有:

{

于是问题得以递归解决: 时间复杂度为 ###### 逆FFT 我们的多项式不可以FFT过去弄不过来了!比较简单, 因为:

所以:

证明略,总之就是证明原向量矩阵可逆所以我们稍加改动就可以逆FFT,就是将单位根的指数取相反数,最后结果再乘1/n 所以两个多项式的乘积,使其对应系数向量的卷积!可以有:

不过我们很少书写递归代码,这里使用迭代写法; 这里还牵扯到两个点:位逆序置换 蝴蝶操作,不说了...直接上代码吧 P3803 【模板】多项式乘法（FFT）

const int N = 3e6+5;
const double pi = acos(-1);
struct Complex
{
    double r, i;
    Complex() {}
    Complex(double r, double i)
    :r(r), i(i) {}
}f[N], g[N];
Complex operator + (const Complex & a, const Complex & b)
{return Complex(a.r + b.r, a.i + b.i);}
Complex operator - (const Complex & a, const Complex & b)
{return Complex(a.r - b.r, a.i - b.i);}
Complex operator * (const Complex & a, const Complex & b)
{return Complex(a.r * b.r - a.i * b.i, a.i * b.r + a.r * b.i);}
int rev[N];
int len, lim = 1;
void init(int n)
{
    while(lim <= n) lim <<= 1, len++;
    for (int i = 0; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len-1));
}
void FFT(Complex *a, int op)//op为1是FFT运算,-1是FFT逆运算
{
    for (int i = 0; i < lim; i++)
    if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int dep = 1; dep <= log2(lim); dep++)
    {
        int m = 1 << dep;
        Complex wn = Complex(cos(2.0 * pi / m), op * sin(2.0 * pi / m));
        for (int k = 0; k < lim; k += m)
        {
            Complex w = Complex(1, 0);
            for (int j = 0; j < m / 2; j++)
            {
                Complex t = w * a[k + j + m / 2];
                Complex u = a[k + j];
                a[k + j] = u + t;
                a[k + j + m / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if (op == -1) for (int i = 0; i < lim; i++)
        a[i].r /= lim;
}
int main()
{
    int n, m;
    cin >> n >> m;
    init(n + m);
    for (int i = 0; i <= n; i++) scanf("%lf", &f[i].r);
    for (int i = 0; i <= m; i++) scanf("%lf", &g[i].r);
    FFT(g, 1); FFT(f, 1);
    for (int i = 0; i <= lim; i++) f[i] = f[i] * g[i];
    FFT(f, -1);
    for (int i = 0; i <= n + m; i++)
        printf("%d ", (int)(f[i].r + 0.5));
    return 0;
}

FFT分治 【模板】分治 FFT 根据cdq分治思想, 我们在计算多项式的一个系数的时候，要统计之前的所有影响。我们区间折半递归解决，先解决左区间的数值，然后统计其对右区间的影响，累加到右区间。洛谷模板需要取模，故模板代码使用NTT，值得注意的是，洛谷模板的模数998244353，符合NTT所要求的素数类型，且存在原根为3。下面是模板代码：

const int N = 1e6+5;
const int G = 3;
const ll mod = 998244353;
ll f[N], g[N];
ll a[N], b[N];
int rev[N], n;
int len, lim = 1;
void init(int n)
{
    lim = 1; len = 0;
    while(lim <= n) lim <<= 1, len++;
    for (int i = 0; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len-1));
}
ll _pow(ll a, ll b) {
    //....
}
inline void calcrev(int logn)
{
    rev[0]=0;
    for(int i = 1; i < (1 << logn); i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(logn-1));
}
void NTT(ll *a, int led, int op)//op为1是FFT运算,-1是FFT逆运算
{
    int lim = 1 << led;
    //......

}
void cdq_fft(int l, int r, int logn)
{
    
    if (logn <= 0) return ;
    if (l >= n) return ;
    int mid = (l + r) >> 1;
    ll t = _pow(r - l, mod-2);
    cdq_fft(l, mid, logn - 1);
    calcrev(logn);
    memset(a + (r - l) / 2, 0, sizeof(ll) * (r - l) / 2);
    memcpy(a, f+l, sizeof(ll) * (r - l) / 2);
    memcpy(b, g, sizeof(ll) * (r - l) );
    NTT(a, logn, 1); NTT(b, logn, 1);
    for (int i = 0; i  < r - l; i++) a[i] = a[i] * b[i] % mod;
    NTT(a, logn, -1);
    for (int i = (r - l) / 2; i < r - l; i++)
        f[l + i] = (f[l+i] + a[i]) % mod;
    cdq_fft(mid, r, logn-1);
}
int main()
{
    cin >> n;
    init(n);
    for (int i = 1; i < n; i++) scanf("%lld", &g[i]);
    f[0] = 1;
    cdq_fft(0, lim, len);
    for (int i = 0; i < n; i++)
        printf("%lld ", (f[i] + mod) % mod);
    puts("");
    return 0;
}