首页 > 其他分享 >多项式学习笔记(一)(2024.7.6)

多项式学习笔记(一)(2024.7.6)

时间:2024-08-01 16:42:24浏览次数:12  
标签:frac 2024.7 int 多项式 sum 笔记 text omega displaystyle

声明:在本节范围内,所有同余号(多项式运算)均在\((\text{mod } x^n)\)意义下进行;所有等号(代数运算)均在模某个质数\(p\)意义下进行。

暴力多项式计算

加法

\(H(x) = F(x) + G(x)\),求\(H(x)\)

解:类比高精度加法

\(h_i = f_i + g_i\),复杂度\(O(n)\)

#include <bits/stdc++.h>
using namespace std;
const int N = 209;
int lenf, f[N];
int leng, g[N];
int lenh, h[N];
int main (){
	scanf("%d", &lenf);
   scanf("%d", &leng);
   lenf++;
   leng++;
	for(int i = lenf; i >= 1; i--)
		scanf("%d", &f[i]);
	for(int i = leng; i >= 1; i--)
		scanf("%d", &g[i]);
	lenh = max(leng, lenf);	
	for(int i = 1; i <= lenh; i++)
		h[i] = f[i] + g[i];   
	for(int i = lenh; i >= 1; i--)
		printf("%d ", h[i]);
	return 0 ;
}

乘法

\(H(x) \equiv F(x)G(x) (\text{mod }x^n)\),求\(H(x)\)

解:类比高精度乘法,复杂度\(O(n^2)\)

#include <bits/stdc++.h>
using namespace std;
const int N = 109;
int f[N], g[N], h[2 * N], lenf, leng, lenh;
int main(){
	scanf("%d", &lenf);
 	scanf("%d", &leng);
 	lenf++;
 	leng++;
	for(int i = lenf; i >= 1; i--)
		scanf("%d", &f[i]);
	for(int i = leng; i >= 1; i--)
		scanf("%d", &g[i]);
	for(int i = 1; i <= lenf; ++i){
        for(int j = 1; j <= leng; ++j){
        	int yx = i + j - 1;
          h[yx] += f[i] * g[j];
        }
	}	
    lenh = lenf + leng;
    for(int i = lenh; i >= 1; --i)
		printf("%d", h[i]);
	return 0;
}

求逆

\(F(x)G(x) \equiv 1 (\text{mod }x^n)\)

求\(G(x)\)

解:设\(H(x) \equiv F(x)G(x)(\text{mod }x^n)\)

\(\therefore h_i = \displaystyle\sum_{j=0}^if_jg_{i-j}\)(多项式乘法)

\(\because H(x)\)除了第\(0\)项是\(1\)其他都为\(0\)。

\(\therefore \displaystyle\sum_{j=0}^if_jg_{i - j} = [i = 0]\)

所有的\(f_i\)已知,那问题就转化成了求一个\(n\)元一次方程组的问题,用高斯消元即可。

除法/取模

\(F(x)\equiv Q(x)G(x)+R(x)(\text{mod }{x^n})\)

求\(Q(x),R(x)\)

解类比高精度除法,复杂度\(O(n^2)\)

#include <bits/stdc++.h>
using namespace std;
const int N = 1000;
int f[N], g[N], q[N], temp[N], lenf, leng, lenq;
int compare(int f[],int g[]){
    if(lenf > leng)
        return 1;
    if(lenf < leng)
        return -1;
    for(int i = lenf; i > 0; i--){  
        if(f[i] > g[i])
            return 1;
        if(f[i] < g[i])
            return -1;
    }
    return 0;
}
void subduction(int f[],int g[]){
    int flag;
    int i;
    flag = compare(f, g);
    if(flag == 0){
        lenf = 0;
        return;
    }
    if(flag == 1){
        for(int i = 1; i <= lenf; i++){
            if(f[i] < g[i]){
                f[i + 1]--;
                f[i] += 10;
            }
            f[i] -= g[i];
        }
        while(lenf > 0 && f[lenf] == 0)
            lenf--;
        return;
    }
}
int main(){
	scanf("%d", &lenf);
	scanf("%d", &leng);
	lenf++;
 	leng++;
	for(int i = lenf; i >= 1; i--)
		scanf("%d", &f[i]);
	for(int i = leng; i >= 1; i--)
		scanf("%d", &g[i]);
    lenq = lenf - leng + 1;
    for(int i = lenq; i > 0; i--){
        memset(temp, 0, sizeof(temp));
        for(int j = 1; j <= leng; j++)
            temp[j + i - 1] = g[j];
        temp[0] = leng + i - 1;
        while(compare(f, temp) >= 0){
            q[i]++;
            subduction(f, temp);
        }
    }
    while(lenq > 0 && q[lenq] == 0)
        lenq--;
    if(lenq == 0)
        printf("0\n");
    else{
        for(int i = lenq; i > 0; i--)
			printf("%d ", q[i]);
        printf("\n");
    }
    if(lenf == 0)
        printf("0\n");
    else{
        for(int i = lenf; i > 0; i--)
			printf("%d ", f[i]);
        printf("\n");
    }
    return 0;  
}

求导/积分

前置知识:求导及其特性

1、定义:\(f'(x) = \mathop{lim}\limits_{\triangle x \rightarrow 0}\displaystyle\frac{f(x + \triangle x) - f(x)}{\triangle x}\)

也可以说导数表示过函数上任意一点的切线斜率

2、性质:①\([f(x) + g(x)]' = \mathop{lim}\limits_{\triangle x \rightarrow 0}\displaystyle\frac{f(x + \triangle x) + g(x + \triangle x)- f(x) - g(x)}{\triangle x} = \mathop{lim}\limits_{\triangle x \rightarrow 0}\displaystyle\frac{f(x + \triangle x) - f(x)}{\triangle x} + \mathop{lim}\limits_{\triangle x \rightarrow 0}\displaystyle\frac{g(x + \triangle x) - g(x)}{\triangle x} = f'(x) + g'(x)\)

②\([f(x)g(x)]'= \mathop{lim}\limits_{\triangle x \rightarrow 0}\displaystyle\frac{f(x + \triangle x)g(x + \triangle x)- f(x)g(x)}{\triangle x} = \mathop{lim}\limits_{\triangle x \rightarrow 0}\displaystyle\frac{f(x + \triangle x)g(x + \triangle x) - f(x + \triangle x)g(x) + f(x + \triangle x)g(x) - f(x)g(x)}{\triangle x} = \mathop{lim}\limits_{\triangle x \rightarrow 0}\displaystyle\frac{f(x + \triangle x)[g(x + \triangle x) - g(x)] + g(x)[f(x + \triangle x) - f(x)]}{\triangle x} = f(x)g'(x) + f'(x)g(x)\)

③复合函数求导:\([f(g(x))]' = f'g(x)g'(x)\)(目前不会证明)

3、特殊函数求导:

①\(f'(x^n) = nx^{n - 1}\)

②\(f'(n^x) = n^x \text{ln }n\)

③\(f'(e^x) = e^x\)

④\(f'(\text{ln } x) = x^{-1}\)

⑤\(f'(\text{sin } x) = \text{cos } x\)

⑥\(f'(\text{cos } x) = -\text{sin } x\)

前置知识:不定积分及其特性

1、定义:若\(g'(x) = f(x)\),则\(\displaystyle\int f(x)dx = g(x) + C\)(加常数是因为求导时会将常数项减掉)

也可以说不定积分表示函数与\(x\)轴之间的面积

2、性质:\(\displaystyle\int [f(x) + g(x)]dx = \displaystyle\int f(x)dx + \displaystyle\int g(x)dx\)(目前不会证明)

3、特殊函数不定积分:

①\(\displaystyle\int f(x^n)dx = \displaystyle\frac{x^{n + 1}}{n + 1} + C\)

②\(\displaystyle\int f(n^x)dx = \displaystyle\frac{1}{\text{ln }a} + C\)

③\(\displaystyle\int f(e^x)dx = e^x + C\)

④\(\displaystyle\int f(\displaystyle\frac{1}{x})dx = \text{ln x} + C\)

⑤\(\displaystyle\int f(\text{sin } x)dx = -\text{cos } x + C\)

⑥\(\displaystyle\int f(\text{cos } x)dx = \text{sin }x + C\)

正题

\(G(x) = F'(x)\)

\(H(x) = \displaystyle\int F(x)dx\)

求\(G(x),H(x)\)

解:\(F(x) = \displaystyle\sum_{i=0}^{n - 1}f_ix^i\)(将多项式展开)

\(F'(x) = \displaystyle\sum_{i=1}^{n - 1}f_i \cdot i \cdot x^{i - 1}\)(特殊函数求导①) \(= \displaystyle\sum_{i=0}^{n - 2}f_{i + 1} \cdot (i + 1) \cdot x^{i}\)(用\(i + 1\)替换\(i\))

又\(\because G(x) = \displaystyle\sum_{i=0}^{n - 1}g_ix^i\)(将多项式展开),第\(x^i\)项系数应该相同

\(\therefore g_i = f_{i + 1} \cdot (i + 1)\)

运用求\(G(x)\)的公式,则有\(f_i = h_{i + 1} \cdot (i + 1)\)

\(\therefore h(i + 1) = \displaystyle\frac{f_i}{i + 1}\)(移项)

\(\therefore h(i) = \displaystyle\frac{f_{i - 1}}{i}\)(用\(i\)替换\(i + 1\))

开根

\(G^2(x) \equiv F(x)(\text{mod }{x^n})\)

求\(G(x)\)

解:运用暴力多项式乘法的公式,则有\(f_i = \displaystyle\sum_{j=0}^{i}g_jg_{i - j} = 2g_0g_i + \displaystyle\sum_{j=1}^{i - 1}g_jg_{i - j}\)

\(\therefore g_i = \displaystyle\frac{f_i - \displaystyle\sum_{j=1}^{i - 1}g_jg_{i - j}}{2g_0}\)

\(\because\)求\(g_i\)时已求出\(g_0\)到\(g_{i - 1}\)

\(\therefore\)可以递推求出\(G(x)\)

注意 \(g_0^2=f_0\) 是边界情况:

  • 若 \(f_0\neq0\),则\(G\)的解数,取决于\(f_0\)的平方根数量。

  • 若 \(f_0=0\),令 \(F(x) = x^kH(x)\),满足 \(h(0) \neq 0\),那么若\(k\)为奇数,则\(G\)无解;若\(k\)为偶数,记 \(P^2(x) = H(x)\),那么 \(G(x) = x^{\frac k2}P(x)\)。

对数函数

\(G(x) \equiv \text{ln }F(x)(\text{mod } x^n)\),保证\(f_0 = 1\)

求\(G(x)\)

解:对原等式两边求导,得:\(G'(x) \equiv F'(x)F^{-1}(x)(\text{mod } x^n)\)(复合函数求导)

令\(H(x) \equiv F^{-1}(x)(\text{mod } x^n)\)(多项式求逆)

\(\therefore G(x) \equiv F'(x)H(x)(\text{mod } x^n)\)

\(\therefore (i + 1)g_{i + 1}=\displaystyle\sum_{j=0}^i{(j + 1)f_{j + 1}h_{i - j}}, g_0 = 0\)

指数函数

\(G(x) \equiv e^{F(x)}(\text{mod } x^n)\),保证\(f_0 = 0\)

求:\(G(x)\)

解:对原等式两边求导,得:\(G'(x) \equiv F'(x)e^{F(x)}\)(复合函数求导) \(\equiv F'(x)G(x)(\text{mod } x^n)\)

\(\therefore (i + 1)g_{i + 1} = \displaystyle\sum_{j=0}^i{(j + 1)f_{j + 1}g_{i - j}}\)

\(\because\)求\(g_{i + 1}\)时已求出\(g_0\)到\(g_i\),且\(g_0 = 1\)

\(\therefore\)可以递推求出\(G(x)\)

三角函数

前置知识:泰勒级数与泰勒展开

\(\because\)求导操作会将常数项消掉,将其它项次数减1。

\(\therefore\)如果两个函数\(F(x), G(x)\)满足\(\left\{ \begin{array}{lr} G(0) = F(0)\\ G'(0) = F'(0)\\ G^{(2)}(0) = F^{(2)}(0)\\ \dots \dots\\ G^{(+\infty)}(0) = F^{(+\infty)}(0) \end{array} \right.\)

那么就可以说明,两个函数相等。(\(G^{(i)}(0) = F^{(i)}(0)\) 实际上表示\(i!g_i = i!f_i\))

\(\therefore\)一个函数\(F(x) = \displaystyle\sum_{i=0}^{+\infty}f_ix^i\)可以表示为\(\displaystyle\sum_{i=0}^{+\infty}\displaystyle\frac{F^{(i)}(0)}{i!}x^i\)

\(\because \text{sin }x = \text{sin }x\)

$\text{sin}'\text{ }x = \text{cos }x $

\(\text{sin}^{(2)}\text{ }x = -\text{sin }x\)

\(\text{sin}^{(3)}\text{ }x = -\text{cos }x\)

\(\text{sin}^{(4)}\text{ }x = \text{sin }x\)(特殊函数求导⑤、⑥)

\(\because \text{sin }x\)在\(x = 0\)处的取值为\(0\)

\(\therefore x^0\)项的系数是\(\displaystyle\frac{0}{0!} = 0\)

\(\because \text{cos }x\)在\(x = 0\)处的取值为\(1\)

\(\therefore x\)项的系数是\(\displaystyle\frac{1}{1!}\)

\(\because \text{-sin }x\)在\(x = 0\)处的取值为\(0\)

\(\therefore x ^ 2\)项的系数是\(\displaystyle\frac{0}{2!} = 0\)

\(\because \text{-cos }x\)在\(x = 0\)处的取值为\(-1\)

\(\therefore x^3\)项的系数是\(\displaystyle\frac{-1}{3!}\)

一直循环下去,则有\(\text{sin }x = \displaystyle\frac 1{1!}x + \displaystyle\frac{-1}{3!}x^3 + \displaystyle\frac1{5!}x^5 + \displaystyle\frac{-1}{7!}x^7 + \dots \dots\)

同理:\(\because \text{cos }x = \text{cos }x\)

$\text{cos}'\text{ }x = -\text{sin }x $

\(\text{cos}^{(2)}\text{ }x = -\text{cos }x\)

\(\text{cos}^{(3)}\text{ }x = \text{sin }x\)

\(\text{cos}^{(4)}\text{ }x = \text{cos }x\)(特殊函数求导⑤、⑥)

\(\because \text{cos }x\)在\(x = 0\)处的取值为\(1\)

\(\therefore x^0\)项的系数是\(\displaystyle\frac{1}{0!}\)

\(\because \text{-sin }x\)在\(x = 0\)处的取值为\(0\)

\(\therefore x\)项的系数是\(\displaystyle\frac{0}{1!} = 0\)

\(\because \text{-cos }x\)在\(x = 0\)处的取值为\(-1\)

\(\therefore x ^ 2\)项的系数是\(\displaystyle\frac{-1}{2!}\)

\(\because \text{sin }x\)在\(x = 0\)处的取值为\(0\)

\(\therefore x^3\)项的系数是\(\displaystyle\frac{0}{3!} = 0\)

一直循环下去,则有\(\text{cos }x = \displaystyle\frac 1{0!}x^0 + \displaystyle\frac{-1}{2!}x^2 + \displaystyle\frac1{4!}x^4 + \displaystyle\frac{-1}{6!}x^6 + \dots \dots\)

同理:\(\because e^{ix} = e^{ix}\)

\((e^{ix})' = ie^{ix}\)

\((e^{ix})^{(2)} = -e^{ix}\)

\((e^{ix})^{(3)} = -ie^{ix}\)

\((e^{ix})^{(4)} = e^{ix}\)(复合函数求导)

\(\because e^{ix}\)在\(x = 0\)处的取值为\(1\)

\(\therefore x^0\)项的系数是\(\displaystyle\frac{1}{0!}\)

\(\because ie^{ix}\)在\(x = 0\)处的取值为\(i\)

\(\therefore x\)项的系数是\(\displaystyle\frac{i}{1!}\)

\(\because -e^{ix}\)在\(x = 0\)处的取值为\(-1\)

\(\therefore x ^ 2\)项的系数是\(\displaystyle\frac{-1}{2!} = 0\)

\(\because -ie^{ix}\)在\(x = 0\)处的取值为\(-i\)

\(\therefore x^3\)项的系数是\(\displaystyle\frac{-i}{3!}\)

一直循环下去,则有\(e^{ix} = \displaystyle\frac 1{0!}x^0 + \displaystyle\frac i{1!}x + \displaystyle\frac{-1}{2!}x^2 + \displaystyle\frac {-i}{3!}x^3 + \displaystyle\frac1{4!}x^4 + \displaystyle\frac i{5!}x^5 + \displaystyle\frac{-1}{6!}x^6 + \displaystyle\frac {-i}{7!}x^7 + \dots \dots = \text{cos }x + i\text{sin }x\)①

于是就有了大名鼎鼎的欧拉公式

于是\(e^{-ix} = \text{cos }x - i\text{sin }x\)②

① \(+\) ② 得:\(2\text{cos }x = e^{ix} + e^{-ix}\)

\(\text{cos }x = \displaystyle\frac{e^{ix} + e^{-ix}}2\)

① \(-\) ② 得:\(2i\text{sin }x = e^{ix} - e^{-ix}\)

\(\text{sin }x = \displaystyle\frac{e^{ix} - e^{-ix}}{2i}\)

正题

\(G(x) = \text{sin }F(x)\)

\(H(x) = \text{cos }F(x)\)

求\(G(x),H(x)\)

直接把\(F(x)\)带进公式,利用多项式指数函数算法即可。

若\(i\)在模\(p\)意义下有对应的数,那么可以直接用这个数来运算。

否则可以扩域,复杂度都是 \(O(n^2)\)。

快速傅里叶变换FFT

前置知识:多项式插值

\(n+1\)个横坐标互不相同的点可以唯一确定一个\(n\)次多项式。

证明:我还不会证,按照两点确定一条直线和三点确定一条抛物线类比吧qwq

前置知识:复数运算

复数是形如 \(a+bi\) 的数,其中 \(i^2=-1\),\(a,b\in\mathbb{R}\)。

复数运算:加:\((a + bi) + (c + di) = (a + c) + (b + d)i\)

乘:\((a + bi)(c + di) = ac - bd + (ad + bc)i\)

除:\(\displaystyle\frac{a + bi}{c + di} = \displaystyle\frac{ac + bd + (bc - ad)i}{c^2 + d^2}\)

前置知识:向量运算

如果将复数\(a + bi\)看成是一条从\((0, 0)\)指向\((a, b)\)的向量,将复数\(c + di\)看成是一条从\((0, 0)\)指向\((c, d)\)的向量,就可以研究向量的加减

向量运算:加:以两条向量为邻边,以两条向量劣角夹角为平行四边形的锐角顶点,做平行四边形,则两向量所夹对角线即为新向量,如图:

乘:将两条向量的模长相乘作为新向量的模长,将两条向量与\(x\)轴的顺时针夹角相加作为新向量与\(x\)轴的顺时针夹角,如图:

前置知识:单位根

声明:以下为了方便书写,均将 \(\omega_n^i\) 表示为 \(\omega_i\)。

代数基本定理:任何复系数一元 \(n\) 次多项式方程在复数域上至少有一根。(我也不会证qwq)

推论:任何复系数一元\(n\)次多项式方程在复数域上恰好\(n\)个根。

证:若一个方程有一根\(x_1\),则此方程可以写成$(x - x_1) \times \((原方程降幂一次的多项式),最多可以降幂\)n\(次,因此恰好有\)n$个根。

单位根:\(x^n - 1 = 0\)的\(n\)个根,记作\(\omega_1, \omega_2, \omega_3 \dots \dots \omega_n\)

求解单位根:在坐标轴中画出单位圆(以\((0, 0)\)为圆心,\(1\)为半径的圆),如图:

将\(x\)当做向量\(\overrightarrow{x}\),则根据向量乘法,\(\overrightarrow{x}\)的模长应为\(1\),由于有\(n\)个根,所以要将这个圆分成\(n\)分,为方便计算,规定\(\overrightarrow{x_1}\)为从\((0, 0)\)指向\((1, 0)\)的向量,就得到了这\(n\)个单位根。以\(n = 6\)时为例子,如图:

由三角函数和复数与向量的关系可得:\(w_i=\cos{\displaystyle\frac{2i\pi}n}+i\sin{\displaystyle\frac{2i\pi}n}\)

如果我们按照这个等式拓宽 \(w_i\) 中 \(i\) 的范围的话,结合其几何意义,我们可以得到 \(w_i\) 一些有趣的性质:

  • 1、\(w_{i+n}=w_i\)(相当于绕着原点转了\(1\)圈)
  • 2、\(w_iw_j=w_{i+j}\)(相当于先顺时针走了\(i\)个单位根,又顺时针走了\(j\)个单位根)
  • 3、\(w_{i}^j=w_{ij}\)(从二性质拓展得来)
  • 4、\(-w_i=w_{-i}\)(顺时针走\(i\)个在取反,等于逆时针走\(i\)个)
  • 5、\(w_{ki}\)=\(w'_{i}\),其中 \(w_i\) 是 \(kn\) 次单位根,\(w'_i\) 是 \(n\) 次单位根。(相当于保持这\(n\)个单位根不变,在每两个单位根之间加入\((k - 1)\)个单位根)

正题:DFT

终于切入正题了,前方高能

根据多项式插值中提到的原理我们可以联想到这样一个解决多项式相乘的算法:

将 \(F(x)\) 和 \(G(x)\) 进行多点求值,得到 \((x_1,F(x_1)),(x_2,F(x_2))\dots(x_n,F(x_n))\) 和 \((x_1,G(x_1)),(x_2,G(x_2))\dots(x_n,G(x_n))\) 这两系列点。

我们可以知道的是:\(H(x)=F(x)G(x)\),所以就有 \(\forall_{i}H(x_i)=F(x_i)G(x_i)\)。

所以我们直接将两系列点对应相乘,得到的一系列点就是 \(H(x)\) 上的点。利用这些点进行插值即可得到 \(H(x)\) 的表达式。

那么如何适当的选取 \(x_1,x_2\dots x_n\),使得多点求值能够快速进行?

我们可以取 \(x_i=w_i\)(\(w_i\)是\(n\)次单位根的第\(i\)个)(因为单位根满足如上的性质)

定义操作\(DFT:\{f_i\}\rightarrow\{F(w_i)\}\),即输入一个\(n-1\)次多项式的系数列,得到其在\(n\)次单位根处的点值列。

现在复杂度还是\(O(n^2)\),需要优化。

当我们想求 \(DFT(F)\) 时,也就是要求 \(F(\omega_k) = \displaystyle\sum_{i=0}^{n - 1}{f_i\omega_k^i}\)

按照奇偶次幂分组,用\(x\)替换\(\omega_k\):原式 \(= (f_0 + f_2x^2 + f_4x^4 + \dots + f_{n - 2}x^{n - 2}) + (f_1x + f_3x^3 + f_5x^5 + \dots + f_{n - 1}x^{n - 1}) = \displaystyle\sum_{i=0}^{\frac n2-1}{f_{2i}x^{2i}} + \displaystyle\sum_{i=0}^{\frac n2-1}{f_{2i+1}x^{2i+1}} = \displaystyle\sum_{i=0}^{\frac n2-1}{f_{2i}x^i} + x\displaystyle\sum_{i=0}^{\frac n2-1}{f_{2i + 1}x^{i}}\)(\(x\)是变量)

设\(F_1(x) = \displaystyle\sum_{i=0}^{\frac n2-1}{f_{2i}x^i}\)。

\(F_2(x) = \displaystyle\sum_{i=0}^{\frac n2-1}{f_{2i + 1}x^{i}}\)

那么可以知道\(F(x) = F_1(x) + xF_2(x)\)

将\(\omega_k(k < \frac n2)\)代入式子中,得到\(F(w_k) = F_1(\omega_k') + \omega_kF_2(\omega_k')\)(\(\omega'\)表示\(n \rightarrow \frac n2\)时的单位根)

将\(\omega_{k + \frac n2}(k < \frac n2)\)代入式子中,得到\(F(w_{k + \frac n2}) = F_1(\omega_{2k}) - \omega_kF_2(\omega_{2k})\)

观察上下两个式子,这两个式子只有一个常数项不同。

那么当我们在枚举第一个式子的时候,我们可以\(O(1)\)的得到第二个式子的值。

容易发现复杂度满足\(T(n)=2T(\frac n2)+O(n)\),根据主定理可知复杂度为 \(O(n\log{n})\)。

于是我们在\(O(n\log n)\)的时间复杂度内求得了 \(DFT(F)\)。

代码:(仅展示\(DFT\))

void DFT(complex <double> *f, int n){
	if(n == 1)
		return;
	for(int i = 0; i < n; i++)//因为是(n - 1)次多项式
		tmp[i] = f[i];
	for(int i = 0; i < n; i++){
		if(i % 2 == 1)
			f[n / 2 + i / 2] = tmp[i];//奇数项放左边
		else
			f[i / 2] = tmp[i];//偶数项放右边
	}
	complex <double> *f1 = f;
	complex <double> *f2 = f + n / 2;
	DFT(f1, n / 2);//递归左边
	DFT(f2, n / 2);//递归右边
	complex <double> cur(1, 0), step(cos(2 * M_PI / n), sin(2 * M_PI / n));//cur记录当前单位根,cur*step让cur成为下一个单位根
	for(int k = 0; k < n / 2; k++){
		tmp[k] = f1[k] + cur * f2[k];//第一个公式
		tmp[k + n / 2] = f1[k] - cur * f2[k];//第二个公式
		cur *= step;//cur变成下一个单位根
	}
	for(int i = 0; i < n; i++)
		f[i] = tmp[i];
}

中置知识:单位根反演

\(\displaystyle\sum_{i=0}^{n - 1}\omega_i^k = n[n \mid k]\)

证:原式 \(= \displaystyle\sum_{i=0}^{n - 1}\omega_{ik}\)(单位根性质3) \(= \displaystyle\sum_{i=0}^{n - 1}\omega_k^i\)(单位根性质3反操作)

现在等式变成了一个等比数列

小知识:等比数列求和

\(q + qp + qp^2 + \dots + qp^n = x\)

解:当\(p = 1\)时,相当于\(n\)个\(q\)相加,\(x = nq\);

当\(p \neq 1\)时,\(\because x = q + qp + qp^2 + \dots + qp^n\)①

\(\therefore px = qp + qp^2 + qp^3 + \dots + qp^{n + 1}\)②

① - ②得:\((1 - p)x = q - qp^{n +1}\)

\(x = \displaystyle\frac{q - qp^{n + 1}}{1 - p}\)


\(\therefore\) 原式 \(= \displaystyle\frac{1 - \omega_k^n}{1 - w_k}\)

当\(n \mid k\)时,有公比 \(w_k = 1\),所以 \(\displaystyle\sum_{i=0}^{n - 1}{\omega_k^i} = n\)

当\(n\nmid k\),有公比 \(w_k \neq 1\),所以 \(\displaystyle\sum_{i=0}^{n - 1}{\omega_k^i} = 0\)

证明完毕。

开始反演:\(nf_k = \displaystyle\sum_{i=0}^{n - 1}w_{-k}^iF(\omega_i)\)

\(F(w_k) = \displaystyle\sum_{i=0}^{n - 1}\omega_k^if_i\)

证:\(nf_k = \displaystyle\sum_{i=0}^{n - 1}\omega_{-k}^i\displaystyle\sum_{j=0}^{n - 1}\omega_k^jf_j\)(将第二个式子代入第一个式子) \(= \displaystyle\sum_{i=0}^{n - 1}\displaystyle\sum_{j=0}^{n - 1}\omega_{-k}^i\omega_k^jf_j\)(将求和移动到最外层) \(= \displaystyle\sum_{i=0}^{n - 1}\displaystyle\sum_{j=0}^{n - 1}\omega_{-ik}\omega_{kj}f_j\)(单位根性质3) \(= \displaystyle\sum_{i=0}^{n - 1}\displaystyle\sum_{j=0}^{n - 1}\omega_{i(j - k)}f_j\)(单位根性质2) \(= \displaystyle\sum_{i=0}^{n - 1}\displaystyle\sum_{j=0}^{n - 1}\omega_i^{j - k}f_j\)(单位根性质3反操作) \(= \displaystyle\sum_{j=0}^{n - 1}f_j\displaystyle\sum_{i=0}^{n - 1}\omega_i^{j - k}\)(交换求和顺序) \(= \displaystyle\sum_{j=0}^{n - 1}f_jn[n \mid (j - k)]\)(第一个公式) \(= \displaystyle\sum_{j=0}^{n - 1}f_jn[j = k]\)(\(j \leq n - 1, k \leq n - 1\),因此\(k - j \leq n\),只有当\(k - j = 0\)时,\((j - k) \div n\)才\(=\)整数,即\(j = k\)) = nf_k

证明完毕

正题:IDFT

定义操作\(IDFT:\{H(w_i)\}\rightarrow\{h_i\}\),即输入一个多项式在\(n\)次单位根处的点值列,得到其\(n-1\)次的系数列。

当我们想求 \(IDFT(h)\) 时,也就是要求 \(nh_k = \displaystyle\sum_{i=0}^{n - 1}{H(\omega_i)\omega_{-k} ^i}\)

有了这个式子之后,我们惊奇的发现\(IDFT\)的操作几乎和\(DFT\)一样,因此同样可以做到\(O(n\log n)\)的复杂度。

代码:(和\(DFT\)很像)

void IDFT(complex <double> *H, int n){
	if(n == 1)
		return;
	for(int i = 0; i < n; i++)
		tmp[i] = H[i];
	for(int i = 0; i < n; i++){
		if(i % 2 == 1)
			H[n / 2 + i / 2] = tmp[i];
		else
			H[i / 2] = tmp[i];
	}
	complex <double> *H1 = H;
	complex <double> *H2 = H + n / 2;
	IDFT(H1, n / 2);
	IDFT(H2, n / 2);
	complex <double> cur(1, 0), step(cos(2 * M_PI / n), sin(2 * M_PI * (-1)/*唯一区别*/ / n));
	for(int k = 0; k < n; k++){
		tmp[k] = H1[k] + cur * H2[k];
		tmp[k + n / 2] = H1[k] - cur * H2[k];
		cur *= step;
	}
	for(int i = 0; i < n / 2; i++)
		H[i] = tmp[i];
}

完整代码:P3803 【模板】多项式乘法(FFT)

#include <bits/stdc++.h>
using namespace std;
const int N = 5e6 + 9;
complex <double> tmp[N << 1];
void DFT(complex <double> *f, int n){
	if(n == 1)
		return;
	for(int i = 0; i < n; i++)
		tmp[i] = f[i];
	for(int i = 0; i < n; i++){
		if(i % 2 == 1)
			f[n / 2 + i / 2] = tmp[i];
		else
			f[i / 2] = tmp[i];
	}
	complex <double> *f1 = f;
	complex <double> *f2 = f + n / 2;
	DFT(f1, n / 2);
	complex <double> cur(1, 0), step(cos(2 * M_PI / n), sin(2 * M_PI / n));
	for(int k = 0; k < n / 2; k++){
		tmp[k] = f1[k] + cur * f2[k];
		tmp[k + n / 2] = f1[k] - cur * f2[k];
		cur *= step;
	}
	for (int i = 0; i < n; i++)
		f[i] = tmp[i];
}
void IDFT(complex <double> *f, int n){
	if(n == 1)
		return;
	for(int i = 0; i < n; i++)
		tmp[i] = f[i];
	for(int i = 0; i < n; i++){
		if(i % 2 == 1)
			f[n / 2 + i / 2] = tmp[i];
		else
			f[i / 2] = tmp[i];
	}
	complex <double> *f1 = f;
	complex <double> *f2 = f + n / 2;
	IDFT(f1, n / 2);
	IDFT(f2, n / 2);
	complex <double> cur(1, 0), step(cos(2 * M_PI / n), sin(2 * M_PI * (-1) / n));
	for(int k = 0; k < n / 2; k++){
		tmp[k] = f1[k] + cur * f2[k];
		tmp[k + n / 2] = f1[k] - cur * f2[k];
		cur *= step;
	}
	for(int i = 0; i < n; i++)
		f[i] = tmp[i];
}
int n, m;
complex <double> f[N << 1], g[N << 1];
int main(){
	scanf("%d%d", &n, &m);
	for(int i = 0; i <= n; i++){
		double tmp;
		scanf("%lf", &tmp);
		complex <double> tmpp(tmp, 0);
		f[i] = tmpp;
	}
	for(int i = 0; i <= m; i++){
		double tmp;
		scanf("%lf", &tmp);
		complex <double> tmpp(tmp, 0);
		g[i] = tmpp;
	}
	int lim = 1;
	while(lim <= n + m)
		lim <<= 1;
	DFT(f, lim);
	DFT(g, lim);
	for(int i = 0; i <= lim; i++)
		f[i] *= g[i];
	IDFT(f, lim);
	for(int i = 0; i <= n + m; i++)
		printf("%d ", int(f[i].real() / lim + 0.5));
	return 0;
}

正题:非递归IDFT

观察IDFT递归顺序,如图:

颜色非常有寓意

最后序列的顺序,等于将原序列下标二进制拆分,前后反转,从新排序。

因此我们可以非递归地完成这一步操作,然后再非递归地往上合并,可以缩减递归所需的花销。

在注意到交换位置的过程实际上可以通过一定的预处理做到\(O(n)\),可以进一步提升效率。

完整代码:P3803 【模板】多项式乘法(FFT)

#include <bits/stdc++.h>
using namespace std;
const int N = 5e6 + 9;
int rev[N];
void DFT(complex <double> *f, int n){
	//找到第一层该数字在什么位置 
	for(int i = 0; i < n; i++){
		rev[i] = rev[i >> 1] >> 1;
		if(i % 2 == 1)
			rev[i] |= n >> 1;
	}
	for(int i = 0; i < n; i++)
		if(i < rev[i])
			swap(f[i], f[rev[i]]);
	for(int h = 2; h <= n; h <<= 1){//每上升一层,每一块数大小加倍 
		complex <double> step(cos(2 * M_PI / h), sin(2 * M_PI / h));
		for(int j = 0; j < n; j += h){//枚举每一块数 
			complex <double> cur(1, 0);
			for(int k = j; k < j + h / 2; k++){//计算H(\omega_k)
				complex <double> c1 = f[k], c2 = cur * f[k + h / 2]; 
				f[k] = c1 + c2;//根据图理解
				f[k + h / 2] =  c1 - c2;
				cur *= step;
			} 
		}
	}
}
void IDFT(complex <double> *f, int n){
	for(int i = 0; i < n; i++){
		rev[i] = rev[i >> 1] >> 1;
		if(i % 2 == 1)
			rev[i] |= n >> 1;
	}
	for(int i = 0; i < n; i++)
		if(i < rev[i])
			swap(f[i], f[rev[i]]);
	for(int h = 2; h <= n; h <<= 1){
		complex <double> step(cos(2 * M_PI / h), sin(2 * M_PI * (-1)/*注意*/ / h));
		for(int j = 0; j < n; j += h){
			complex <double> cur(1, 0);
			for(int k = j; k < j + h / 2; k++){ 
				complex <double> c1 = f[k], c2 = cur * f[k + h / 2]; 
				f[k] = c1 + c2;
				f[k + h / 2] =  c1 - c2;
				cur *= step;
			} 
		}
	}
}
int n, m;
complex <double> f[N << 1], g[N << 1];
int main(){
	scanf("%d%d", &n, &m);
	for(int i = 0; i <= n; i++){
		double tmp;
		scanf("%lf", &tmp);
		complex <double> tmpp(tmp, 0);
		f[i] = tmpp;
	}
	for(int i = 0; i <= m; i++){
		double tmp;
		scanf("%lf", &tmp);
		complex <double> tmpp(tmp, 0);
		g[i] = tmpp;
	}
	int lim = 1;
	while(lim <= n + m)
		lim <<= 1;
	DFT(f, lim);
	DFT(g, lim);
	for(int i = 0; i <= lim; i++)
		f[i] *= g[i];
	IDFT(f, lim); 
	for(int i = 0; i <= n + m; i++)
		printf("%d ", int(f[i].real() / lim + 0.5));
	return 0;
}

快速数论变换NTT

前置知识:原根

若\(a^{\varphi(p)} \equiv 1 (\text{mod } m)\), 且对\(\forall x \in [0, \varphi(p) - 1], a^x \not\equiv 1 (\text{mod } m)\),则称\(a\)为\(p\)的原根。

有趣的知识:\(998244353\)的一个原根是\(114514\)(°O °)

正题

设\(g\)表示\(p\)的一个原根,则令 \(\omega_i=g^{i\frac{\varphi(p)}n}\) 就具有同复数域中单位根一样的性质。(见快速傅里叶变换FFT/前置知识:单位根/单位根的性质)

但是这不代表我们可以直接将他套用进 FFT 中。因为这里必须要求 \(n\mid \varphi(p)\),所以这对 \(n\) 产生了很大的限制。

幸运的是我们发现质数 \(998244353\) 的欧拉函数值 \(998244352=119\times2^{23}\),因此当\(n \leq 2^{23}\) 时,都能够很好的进行这个算法,这也是\(998244353\)成为常用模数的原因。(尽管有些时候题目与多项式无关)

完整代码:P3803 【模板】多项式乘法(FFT)

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 5e6 + 9, MOD = 998244353;
int rev[N];
int qpow(int x, int y){
	int res = 1;
	while(y){
		if(y & 1)
			res = res * x % MOD;
		x = x * x % MOD;
		y >>= 1;
	}
	return res;
}
void NTT(int *f, int n, int opt){
	for(int i = 0; i < n; i++){
		rev[i] = rev[i >> 1] >> 1;
		if(i % 2 == 1)
			rev[i] |= n >> 1;
	}
	for(int i = 0; i < n; i++)
		if(i < rev[i])
			swap(f[i], f[rev[i]]);
	for(int h = 2; h <= n; h <<= 1){
		int mi = h >> 1;
		int gn = qpow(3, (MOD - 1) / h);
		for(int j = 0; j < n; j += h){
			int g = 1;
			for(int k = 0; k < mi; k++){
				int tmp = f[j + k + mi] * g % MOD;
				f[j + k + mi] = (f[j + k] - tmp + MOD) % MOD;
				f[j + k] = (f[j + k] + tmp) % MOD;
				g = g * gn % MOD;
			}
		}
	}
	if(opt == -1){
		reverse(f + 1, f + n);
	    int inv = qpow(n, MOD - 2);
	    for(int i = 0; i < n; i++)
			f[i] = f[i] * inv % MOD;
	}
}
int f[N << 1], g[N << 1], n, m;
signed main(){
	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 lim = 1;
	while(lim <= n + m)
		lim <<= 1;
	NTT(f, lim, 1);
	NTT(g, lim, 1);
	for(int i = 0; i <= lim; i++)
		f[i] = f[i] * g[i] % MOD;
	NTT(f, lim, -1); 
	for(int i = 0; i <= n + m; i++)
		printf("%lld ", f[i]);
	return 0;
}

标签:frac,2024.7,int,多项式,sum,笔记,text,omega,displaystyle
From: https://www.cnblogs.com/JPGOJCZX/p/18336949

相关文章

  • [学习笔记]最小割树 (Gomory-Hu Tree)
    本文是在作者阅读多篇博客后糅合而成的,部分内容有雷同。最小割树又称Gomory-Hu树,由RalphEdwardGomory和TeChiangHu于1961年发表的论文中共同提出。最小割树可以在\(n−1\)次最大流中构建,并通过树上RMQ求任意两点之间的最小割。板子题:P4897【模板】最小割树(G......
  • ScottPlot使用笔记
    安装WPF:Install-PackageScottPlot.WPFAvalonia:Install-PackageScottPlot.Avalonia相关网站https://github.com/ScottPlot/ScottPlothttps://scottplot.net/调研情况ScottPlot调研情况:目前可以用ScottPlot来绘制基本的曲线图,柱形图,饼图,雷达图。ScottPlot最新的是5.X......
  • MySQL 学习笔记 进阶(InnoDB引擎 下)
    InnoDB引擎 InnoDB引擎-事务原理-概述事务是一组操作的集合,它是一个不可分割的工作单位,事务会把所有的操作作为一个整体一起向系统提交或撤销操作请求,即这些操作要么同时成功,要么同时失败。原子性(Atomicity):事务是不可分割的最小操作单元,要么全部成功,要么全部失败。一致性(Co......
  • 【笔记】字符串选讲 2024.8.1
    [COCI2015-2016#5]OOP(Trie)P6727[COCI2015-2016#5]OOP-洛谷|计算机科学教育新生态(luogu.com.cn)正反串分别建Trie,可以搞出两个dfn区间,加之长度限制,三维数点。有\(O(n\logn)\)做法。将字典串\(S[1..m]\),对所有\(1\leqi\leqm\),将\(S[i+1,m]\)的hash值插入......
  • Linux操作系统基础学习笔记(4)
    Linux操作系统基础学习笔记(4)前言4、Linux文件和目录管理常规命令格式(1)列出目录内容和属性(文件)(2)打印工作路径(3)切换工作路径(4)查看文件类型(5)复制文件或目录(6)查找文件或目录(7)创建目录(8)移动或重命名(9)删除文件(不能用来删除文件夹)(10)创建空文件(11)挂载(12)链接(有点像windows的快捷......
  • 组合数学学习笔记(持续完善中)
    基础知识一、加法原理完成某个工作有\(n\)类办法,第\(i\)类办法有\(a_i\)种,则完成此工作的方案数有\(\sum\limits_{i=1}^na_i\)种。二、乘法原理完成某个工作有\(n\)个步骤,第\(i\)个步骤有\(b_i\)种,则完成此工作的方案数有\(\prod\limits_{i=1}^nb_i\)种。......
  • 努力努力努力的第十四天(2024.7.31)
    昨天日期写错了写成2020.7.30,应该是2024.7.31(手滑了哈哈哈)1.行列转换效果演示:这是未经行列转换操作的t_score表:这是经过行列转换后的t_score表:第一步:确定初步的做法使用分组查询(groupby)能够将单个学生的成绩依次查询出来,再加上三列查询(分别定义成'语文''数学''......
  • java笔记2
    4.数组笔记数组概念数组是一种基本的数据结构,用于存储固定大小的相同类型的元素序列。在Java中,数组是一种对象,它实现了java.lang.Cloneable和java.io.Serializable接口。声明数组:int[]intArray;初始化数组:intArray=newint[10];//创建一个长度为10的整型数组......
  • 【远程驰骋:Python SSH 自动化运维实战笔记】
    使用GqylpySSH库简化SSH命令执行在自动化运维或脚本编写中,经常需要通过SSH连接到远程服务器执行命令。虽然Python的paramiko库提供了强大的SSH功能,但直接使用它进行命令执行和结果处理可能会显得有些繁琐。GqylpySSH库封装了paramiko,提供了一个更加简洁易用的接口......
  • FreeRTOS学习笔记(二)
    FreeRTOS移植一、获取FreeRTOS源码1.1官网下载进入官网直接下载官网:https://www.freertos.org/zh-cn-cmn-s/1.2正点原子网盘下载正点原子资料v10.4.6例程git:https://gitee.com/yuan-zhenbin/freertos-code-repository.gitFreeRTOS资料网盘:http://www.openedv.c......