首页 > 其他分享 >【多项式】【模版】FFT NTT

【多项式】【模版】FFT NTT

时间:2024-02-16 20:22:42浏览次数:22  
标签:dots int 模版 FFT NTT 复数 多项式 点值

多项式

若 \(A(x)=a_0+a_1x+a_2x^2+ \dots a_nx^n(a_n \ne 0)\),则称 \(A\) 是 \(n\) 次多项式。

初中概念,但在 OI 中可以玩出花来。

多项式的表示方式

像上面的介绍一样,用系数 \(a_0,a_1, \dots a_n\) 来表示多项式的方法称为系数表示法

还有一种表示多项式的方法,就是对于 \(n\) 次多项式,给出 \(n+1\) 对点,\((x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2)) \dots (x_n,A(x_n))\), 且 \(x_0,x_1 \dots x_n\) 互异,称为 点值表示法

可以证明,点值表示法和系数表示法一一对应,表示唯一的多项式。

系数表示法 \(\rightarrow\) 点值表示法:带入计算即可,复杂度 \(O(n^2)\)。

点值表示法 \(\rightarrow\) 系数表示法:需要使用拉格朗日插值,复杂度 \(O(n^2)\)。

之后我们引入 FFT,可以实现 \(O(n \log n)\) 时间内两种表示法的转化。

多项式乘法

多项式乘法是多项式的核心操作。


若有两个多项式 \(A(x) = \sum_{i=0}^na_ix^i\),\(B(x) = \sum_{i=0}^mb_ix^i\),其乘积为

\[M(x)=A(x)B(x)=\sum_{i=0}^{n}\sum_{j=0}^{m}a_ib_jx^{i+j} \]

就是简单的代数运算,\(n\) 次多项式和 \(m\) 次多项式的乘积次数是 \(n+m\)。

直接硬乘复杂度自然是 \(O(n^2)\)。


对于点值表示法,两个多项式预先准备 \(n+m+1\) 对点,然后乘起来,得 \((x_{ai}x_{bi},A(x_{ai})B(x_{b_i}))\)。

这 \(n+m+1\) 对点足够表示出多项式。 复杂度是 \(O(n)\)。

可以看出,点值表示法做乘法速度很快,但系数表示法无论是直接乘还是转化为点值乘,都是 \(O(n^2)\),为了解决这个问题,就要请出 FFT——快速傅里叶变换了

FFT

简介

我们如果将 \(n+1\) 个普通的数带入多项式计算点值,复杂度难以优化。而 FFT,简单来说,就是通过带入特殊的数,来加速系数表示 \(\rightarrow\) 点值表示的速度

单位复数根

基础知识:数学必修二(其实就是复数)

在复数域内,对于方程 \(w^n=1\),一共有 \(n\) 个解,由

\[e^{iu}=\cos u+i \sin u \]

可得,为 \(e^{2 \pi ik/n}\),其中 \(k=0,1\dots n-1\),称为 \(n\) 次单位复数根。

对于 \(k=1\),我们定义其为 \(w_n=e^{2\pi i/n}\),称为主 \(n\) 次单位复数根。那么这 \(n\) 个根就可以表示为 \(w_n^0,w_n^1 \dots w_n^{n-1}\),它们均匀分布在复平面单位圆上。

性质一

\(w_{dn}^{dk}=w_{n}^{k}\)(因为单位复数根是在单位圆上均分,所以只看比例 \(k/n\))

引理

若 \(n\)为偶数,\((w_n^k)^2=w^k_{n/2}\)

性质二

若 \(n\) 为偶数,\(w_n^k=-w_n^{k+n/2}\)(几何解释:转半圈)

快速傅里叶变换

利用单位复数根的性质,我们通过分治法代入 \(n\) 次单位复数根,来得到 \(n-1\) 次多项式的点值表示。

基于这样一个分治思想,令 \(A^{[0]}(x) = a_0+a_2x+a_4x^2\dots\),\(A^{[1]}(x)=a_1+a_3x+a_5x^2\dots\),就是把偶数项系数和奇数项系数分给两个多项式,那么就有

\[A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2) \]

对 \(A^{[0]}(x),A^{[0]}(x)\) 分治计算即可。

这样并不能优化复杂度,因为普通的值代入多项式,\(A^{[0]},A^{[1]}\) 实际要计算 \(x_1, x_2^2,x_3^2\dots\) 共 \(n\) 个值,但又有它们是 \(n/2\) 次多项式,理论上只需要代入 \(n/2\) 个值。

这个时候,单位复数根就起作用了,如果计算 \(A(w_n^k)\),分治时,\((w_n^k)^2\) 变为了 \(w_{n/2}^k\),就有

\[A(w_n^k)=A^{[0]}(w_{n/2}^k)+w_n^kA^{[1]}(w_{n/2}^k) \]

同时根据性质二,有

\[A(w_n^{k+n/2})=A^{[0]}(w_{n/2}^k)-w_n^kA^{[1]}(w_{n/2}^k) \]

那么对于 \(A^{[0]},A^{[1]}\),它们只要计算 \(n/2\) 个值就行了,复杂度为 \(T(n)=2T(n/2)+n\),为 \(O(n \log n)\)。

边界是 \(n=1\) 的情况,此时 \(A(w_1^0)=A(1)=a_0\)

下面是 FFT 的递归模版,将系数表示原址转化为点值表示

void FTT(complex<double>* A, int n) {
    if (n == 1) return ; // A(w_1^0)=A(1)=a_0
    vector<complex<double>> tmp(n);
    complex<double> wn(cos(2 * pi / n), sin(2 * pi / n)); //n次单位复数根计算
    complex<double> wk(1, 0);
    for (int i = 0; i < n; i++) {
        if (i & 1) tmp[i / 2 + n / 2] = A[i];
        else tmp[i / 2] = A[i];
    }
    for (int i = 0; i < n; i++) A[i] = tmp[i]; 
    //  对A数组进行原址划分,递归左右两部分
    // 得出来的A数组[0,n/2)内储存A0(w_n/2)
    //            [n/2,n)内储存A1(w_n/2)
    vector<complex<double>>().swap(tmp);
    FTT(A, n / 2), FTT(A, n / 2);
    for (int i = 0; i < n / 2; i++, wk *= wn) {
        auto x = A[i], y = A[i + n / 2]; 
        A[i] = x + wk * y;
        A[i + n / 2] = x - wk * y;
    }
}

在 FFT 过程中,要确保 \(n\) 随时为偶数,所以最初 \(n\) 必须得是 \(2\) 的幂。

这里使用了 c++ 的复数模版进行运算,如果只在理论层面上学习 FFT,已经可以完结撒花了。但是,这个模版有几个缺陷。

1.使用复数运算,递归,且多次复制数组,造成了时间空间的巨大开销/

解决方案:自定义复数运算,学习接下来的FFT的迭代实现。

2.三角函数带来的精度问题

解决方案:我们之后介绍的快速数论变换会解决整数多项式在精度上的问题。

FFT 的迭代实现

为了优化常数,我们要将递归改为循环迭代,这是个很考验 FFT 原理的活。

首先递归FFT中有“划分奇偶系数”这一操作,观察递归树。

原来的 \(a_0,a_1\dots a_8\) 变为了 \(a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7\),如果细心的话,可以发现最后 \(a_i\) 所在的下标就是
\(i\) 二进制倒过来,举两个例子

\(a_3:011 \rightarrow 110; 3\rightarrow 7\)

\(a_4:100 \rightarrow 001; 4\rightarrow 1\)

证明很简单,第 \(k\) 层划分是看 \(i\) 第 \(k\) 位是否为 \(1\),然后以 \(2^{dep-k}\) 的权值加到新下标上,所以就是二进制倒过来啦。

我们可以 \(O(n)\) 递推求解 \(0\dots n-1\) 内的倒二进制。

for (int i = 0; i < n; i++) { //自行模拟
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1) rev[i] |= n >> 1;
}

这被称为 位逆序置换


“划分”完了,我们要考虑“合并”,这个和递归的合并并无二异,当前 \(A^{[0]}(x)\) 和 \(A^{[1]}(x)\) 分别记录在 \([0,n/2)\) 和 \([n/2,n)\) 中。调用然后覆盖即可。

下面给出代码

//已经计算好 rev
void FFT(auto *A, int n) {
    for (int i = 0; i < n; i++) {
        if (i < rev[i]) swap(A[i], A[rev[i]]);
    }
    for (int i = 1; i < n; i <<= 1) { //这里i是枚举 n/2,
        auto wn = ...
        for (int j = 0; j < len; j += (i << 1)) { //处理[j,j+n)
            auto wk(1, 0)...
            for (int k = 0; k < i; k++, wk = wk * wn ) {
                auto x = A[j + k], y =  A[i + j + k];
                A[j + k] = x + wk * y;
                A[i + j + k] = x - wk * y;
            }
        }
    }
}

逆变换——IDFT

我们由上文知道了系数转点值,那么点值怎么转系数呢?

注意,这里的点值应是通过 FFT 转换过来的 \(A(w_n^k)\),

直接看结论:对点值数组以 \(w_n^{-k}\) 为单位复数根,进行一次 FFT,然后每一个数除以 \(n\)。

我们直接把单位复数根判断写进 FFT 里,就能 FFT IDFT 一起了。

证明:

有 \(A(w_n^k)=a_0+a_1w_n^k+a_2w_n^{2k}+\dots\)

\(A(w_n^k)w_n^{-ki}=a_0w_n^{-ki}+a_1w_n^{k(1-i)}+a_2w_n^{k(2-i)}+\dots\)

我们对 \(k\) 求和:

\(\displaystyle\sum_{k=0}^{n-1}{A(w_n^k)w_n^{-ki}}=a_0(w_n^0+w_n^{-i}+w_n^{-2i}+\dots)+a_1(w_n^0+w_n^{1-i}+w_n^{2(1-i)}+\dots)+\dots\)

由等比数列求和公式,可得:

\(w_n^{0k}+w_n^{1k}+\dots w_n^{(n-1)k}=\displaystyle\frac{w_n^{nk}-w_n^0}{w_n^k-1}=0\)

所以对于第 \(i\) 项 \(w_n^{-i}\),将其代入多项式 \(F(x)=A(w_n^0)+A(w_n^1)x+A(w_n^2)x^2+\dots\) 中,非 \(i\) 项全部被消掉,只留下 \(na_i\),所以就是这样了。

快速数论变换——NTT

这个其实可以单讲的,但其实原理和 FFT 一模一样,只不过把单位复数根换成了原根。对于整数多项式,消除了误差。

前置知识:原根(这玩意真的很有用啊)

我们使用单位复数根是因为它很好的性质,平方后数量减半,我们用原根也能达到这样的效果。由于简便,下面的等号都是在模意义上相等。

选取一个奇素数 \(p\),然后找到原根 \(g\),

我们令 \(n\) 次主单位根 \(w_n=g^{(p - 1)/n}\),可以发现 \(w_n^n=g^{p-1}=1=w_n^{0}\),形成了一个环。

然后对于性质一,自然成立;

对于性质二,\(w_n^{n/2}= g^{(p-1)/2}=-1\),所以成立。

那么和 FFT 代码基本一样,只不过把复数换成了模运算。

由于 \(n\) 要整除 \(p-1\),还得是 \(2\) 的幂,我们尽可能让 \(p-1\) 有较多的 \(2\),可以选取 \(998244353=2^{23}*119+1\)。

最终模版(我还是觉得NTT好用)

const int N = 1 << 22;
const int mod = 998244353, g = 3;
const int gi = pow_mod(3, mod - 2);
void init() {
    for (int i = 0; i < len; i++) {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1) rev[i] |= len >> 1;
    }
}
void NTT(int *A, int len, bool t) {
    for (int i = 0; i < len; i++) {
        if (i < rev[i]) swap(A[i], A[rev[i]]);
    }
    for (int i = 1; i < len; i <<= 1) {
        int wn = pow_mod(t ? g : gi, (mod - 1) / (i << 1));
        for (int j = 0; j < len; j += (i << 1)) {
            int w0 = 1;
            for (int k = 0; k < i; k++, w0 = w0 * wn % mod) {
                int x = A[j + k], y = w0 * A[i + j + k] % mod;
                A[j + k] = (x + y) % mod;
                A[i + j + k] = (x - y + mod) % mod;
            }
        }
    }
    if (!t) {
        int inv = pow_mod(len, mod - 2);
        for (int i = 0; i < len; i++) A[i] = A[i] * inv % mod;
    }
}

标签:dots,int,模版,FFT,NTT,复数,多项式,点值
From: https://www.cnblogs.com/Uuuuuur/p/18017442

相关文章

  • Django 模版变量
    Django模版变量Django模版语言的语法主要分为以下四个部分:变量标签过滤器注释一、模版变量1)变量的命名规范Django对于模版变量的命名规范没有太多的要求,可以使用任何字母、数字和下划线的组合来命名,且必须以字母或下划线开头,但是变量名称中不能有空格或者标点符号。2......
  • 十二、Django视图函数和模版相关
    视图相关HTTPRequest对象:属性:path函数:get_full_path()HTTPResponse对象:render()render_to_response()locals():局部变量redirect()重定向例子:用户登录成功后跳转deflogin(request):...#判断登录成功后,跳转到indexreturnrender(request,"index.html"......
  • P4721 【模板】分治 FFT
    最具经济性的写法:\(\mathcalO(n^2)\)暴力拿下\(80\)分,遂跑路。一题多解了,分两部分:分治和多项式求逆。分治考虑cdq分治,每次把\(f_{l\dotsmid}\)和\(g_{1\dotsn-1}\)卷起来,贡献直接加到\(f_{mid+1\dotsr}\)里,要注意一下顺序,先递归左区间,再算当前区间,最......
  • 【数论】【模版】二次剩余
    二次剩余问题其实就是数论中的开方运算。我们要解决这么一个问题,给定正整数\(n\),奇素数\(p\),求解\[x^2\equivn\pmodp\]本文内认为模数\(p\)是一个奇素数。定义若存在\(x^2\equivn\pmodp\),则称\(n\)为模\(p\)的二次剩余,反之则称\(n\)为模\(p\)的非二次剩......
  • 3秒钟教你如何配置vscode中的vue3代码快速生成模版
    1.首先点击你的vscode左下角的齿轮设置按钮,然后点击配置用户代码片段2.输入vue搜索vue.json这个文件,然后点击这个文件3.接下来只需在原有的注释之下输入粘贴如下代码即可4.代码如下"vue3":{"prefix":"vue3","body":["<template>",......
  • FFT快速傅里叶变换
    傅里叶变换复数定义:x,y为实数,称形如(x,y)的有序数对为复数。复数(x,y)中的第一个实数x称为复数z的实部,第二个实数y称为复数z的虚部。代码实现及运算:structComplex{ doux,y; Complex(douxx=0,douyy=0){x=xx,y=yy;} Complexoperator+(Com......
  • FFT理论与习题
    参考了这篇博客,并且自己重新证了一下这篇博客中,笔者认为有误的IDFT中\(j\neqk\)的部分。第零部分·原理FFT是一种DFT的高效算法,称为快速傅里叶变换(fastFouriertransform),当然也可以用来算IDFT。这俩玩意前者可以把多项式从系数表达式转化成点值表达式,后者可以把点......
  • 上辈子推的"分治 NTT"复杂度分析
    hdu7381Cargo式子部分由liuhangxin想出。\[\sum\limits_{i=0}^{n}\binom{k}{i}(n-i)^{k-i}[x^i]\prod\limits_{id=1}^{m}(1-x^{c_{id}})\]实现部分当时胡了一个分治NTT,也不知道时间复杂度为什么是对的,但是过了。AC后十多分钟分析出来这个做法的时间复杂度为\(......
  • 【快速阅读四】基于边缘信息的模版匹配中贪婪度参数的简单解析。
    基于边缘的匹配有个贪婪度的参数,其对提高查找目标的速度有着比较关键的作用,本问简单的记录了下对这个参数的一些认识和推导。对这个课题稍作研究,以便记录。在基于边缘的模版匹配中,我们知道可以有个贪婪度参数可以设置。在Halcon的帮助文档中,也有对......
  • audio currentTime 设置后,重置成0,解决方案(流文件-下载文件)
    audiocurrentTime设置后,重置成0,解决方案第一步-流文件-下载文件:先查看你的mp3文件是流文件,还是下载文件。检测方式,就是放到浏览器回车。在线播放就是流文件,直接下载了,就是下载文件。是需要流文件,才可以拖拽进度条。第二步:currentTimeseekTo方法针对不同浏览器,做下兼容......