工科数学分析大作业(三) 傅里叶级数
研究课题:探索类问题7
给定多项式:(1)∑n=023n2xn,∑n=024nxn (1) \sum\limits_{n=0}^{2^3} n^2 x^n , \sum\limits_{n=0}^{2^4} n x^n ; (2)∑n=0220n2xn,∑n=0220nxn (2) \sum\limits_{n=0}^{2^{20}} n^2 x^n , \sum\limits_{n=0}^{2^{20}} n x^n ;
利用快速傅里叶变换实现多项式相乘,并分析算法的复杂度。
算法介绍:快速数论变换
模 PP 意义下的多项式
A(x)=∑i=0n−1aixiA(x)=\sum\limits_{i=0}^{n-1} a_i x^i ,称 a0a_0 , a1a_1 , ⋯\cdots , an−1a_{n-1} 为多项式的系数。所有系数都属于模 PP 的完全剩余系。如果多项式 A(x)A(x) 的最高次的非零系数为 aka_k ,则称 A(x)A(x) 的次数是 kk . 任何一个严格大于一个多项式次数的整数都是这个多项式的次数界。
多项式乘法
如果多项式 A(x)A(x) 和 B(x)B(x) 都是次数界为 nn 的多项式,则说它们的乘积是一个次数界为 2n−12n-1 的多项式积 C(x)C(x),并满足对所有属于定义域的 xx,都有 C(x)=A(x)⋅B(x)modPC(x)=A(x) \cdot B(x) \; mod \; P .
多项式的表示方法
系数表示法:
对一个次数界为 nn 的多项式 A(x)=∑i=0n−1aixiA(x)=\sum\limits_{i=0}^{n-1} a_i x^i 来说,其系数表示法就是一个由系数组成的向量 a=(a0,a1,⋯,an−1)(a_0,a_1,\cdots,a_{n-1}) 。
系数表示法下的多项式乘法
假设 AA 和 BB 都是次数界为 nn 的多项式,aa 和 bb 分别是 AA 和 BB 的系数向量,则 AA 和 BB 的乘积 CC 的系数向量 c=a⨂b=(c0,c1,⋯.cn−1)c=a\bigotimes b=(c_0,c_1,\cdots.c_{n-1}),其中 ci≡∑j=0iaj⋅bi−jmodPc_i\equiv\sum\limits_{j=0}^{i} a_j \cdot b_{i-j}\;mod\;P 。根据这个定义计算两个次数界为 nn 的多项式乘法的时间复杂度为 O(n2)O(n^2) 。
点值表示法:
一个次数界为 nn 的多项式 A(x)A(x) 的点值表示就是 nn 个点值对所形成的集合:{(x0,y0),(x1,y1),⋯,(xn−1,yn−1)}\{ (x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1}) \},其中所有 xkx_k 互不相同,并且当 k=0,1,⋯,n−1k=0,1,\cdots,n-1 时有yk=A(xk)modPy_k= A(x_k)\;mod\;P。一个多项式可以有很多不同的点值表示。
如果任意选取 {xk}\{x_k\},根据霍纳法则(秦九韶算法)将 A(x)A(x) 从系数表示法转换为点值表示法(求值)的时间复杂度为 O(n2)O(n^2) 。如果巧妙的选取 {xk}\{x_k\},可以将转换的时间复杂度降为 O(nlogn)O(nlogn) 。
求值计算的逆(从一个多项式的点值表示确定其系数表示中的系数)称为插值。
多项式插值的唯一性
对于任意 nn 个相异点值对组成的集合 {(x0,y0),(x1,y1),⋯,(xn−1,yn−1)}\{(x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\},存在唯一的次数界为 nn 的多项式 A(x)A(x),满足 yk≡A(xk)modPy_k\equiv A(x_k)\;mod\;P,k=0,1,⋯,n−1k=0,1,\cdots,n-1 。
证明: 证明过程的基础是某个矩阵存在逆矩阵,yk=A(xk)modPy_k=A(x_k)\;mod\;P,k=0,1,⋯,n−1k=0,1,\cdots,n-1 等价于矩阵方程
⎡⎣⎢⎢⎢⎢⎢11⋮1x0x1⋮xn−1x20x21⋮x2n−1⋯⋯⋱⋯xn−10xn−11⋮xn−1n−1⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢y0y1⋮yn−1⎤⎦⎥⎥⎥⎥⎥\begin{equation} \left[ \begin{matrix}1&x_0&x_0^2&\cdots&x_0^{n-1}\\1&x_1&x_1^2&\cdots&x_1^{n-1}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&x_{n-1}&x_{n-1}^2&\cdots&x_{n-1}^{n-1} \end{matrix} \right] \left[ \begin{matrix} a_0\\a_1\\\vdots\\a_{n-1} \end{matrix} \right] = \left[ \begin{matrix} y_0\\y_1\\\vdots\\y_{n-1} \end{matrix} \right] \end{equation}
左边的矩阵是范德蒙德矩阵,表示为 V(x0,x1,⋯,xn−1)V(x_0,x_1,\cdots,x_{n-1})。该矩阵的行列式的值为 ∏0≤i<j≤n−1(xj−xi)\prod\limits_{0\leq i,因此,如果 xkx_k 相异,则该矩阵是可逆的。因此,对给定的唯一的点值表示,就能够求出系数 aia_i 的解: a=V(x0,x1,⋯,xn−1)−1ya=V(x_0,x_1,\cdots,x_{n-1})^{-1}y 。
点值表示法下的多项式乘法
点值表示法对于多项式乘法是很方便的。如果 C(x)=A(x)B(x)C(x)=A(x)B(x),则对任意点 xkx_k,有 C(xk)=A(xk)B(xk)C(x_k)=A(x_k)B(x_k),并且对 AA 的点值表示和 BB 的点值表示进行逐点相乘,就可以得到 CC 的点值表示。假设 AA 和 BB 的次数界都为 nn,则 CC 的次数界为 2n2n,因此要得到 CC 的点值表示需要 2n2n 个点值对。因此,必须对 AA 和 BB 的点值表示进行“扩充”,使每个多项式都包含 2n2n 个点值对。
如果已知 AA 的扩充点值表示:{(x0,y0),(x1,y1),⋯,(x2n−1,y2n−1)}\{(x_0,y_0),(x_1,y_1),\cdots,(x_{2n-1},y_{2n-1})\} 和 BB 的相应扩充点值表示:{(x0,y′0),(x1,y′1),⋯,(x2n−1,y′2n−1)}\{(x_0,y^{'}_0),(x_1,y^{'}_1),\cdots,(x_{2n-1},y^{'}_{2n-1})\},则 CC 的点值表示为:{(x0,y0y′0modP),(x1,y1y′1modP),⋯,(x2n−1,y2n−1y′2n−1modP)}\{(x_0,y_0y^{'}_0\;mod\;P),(x_1,y_1y^{'}_1\;mod\;P),\cdots,(x_{2n-1},y_{2n-1}y^{'}_{2n-1}\;mod\;P)\} 。
如果已知两个扩充点值形式的输入多项式,使其相乘而得到点值形式的结果需要 O(n)O(n) 的时间。
原根
设 mm 是正整数,aa 是整数,若 aa 模 mm 的阶等于 φ(m)\varphi(m) ,则称 aa 为模 mm 的一个原根(其中 φ(m)\varphi(m) 是 mm 的欧拉函数)。
假设 gg 是 PP 的一个原根,则 g0,g,g2,⋯,gP−2g^0,g,g^2,\cdots,g^{P-2}在模 PP 的乘法运算下形成了一个 P−1P-1 阶循环群,其中 PP 为素数,且有 gP−1≡g0≡1modPg^{P-1} \equiv g^0 \equiv 1 \;mod\;P 。
系数形式表示的多项式的快速乘法
能否采用利用关于点值形式的表示的多项式的线性时间乘法算法,来加快系数形式表示的多项式乘法运算的速度呢?答案依赖于能否快速把一个多项式在系数形式和点值形式之间转换。
当模数特殊的时候,当 P=r2k+1P=r2^k+1,且 PP 为素数,存在原根 gg 的时候,如果选取 gP−1ng^{\frac{P-1}{n}} 作为单位根,则可以把一个次数界为 nn 的多项式在两种表示法之间转化所需的时间压缩为 O(nlogn)O(nlogn) (这里要求 n=2jn=2^j 并且 j≤kj\leq k)。
本文中之后所提到的所有模数 PP 均默认为满足上述条件。
NTT和FNTT
模P意义下的单位根
类似于复数域上单位复根的定义,在模 PP 的意义下定义 nn 次单位根为满足 ωn≡1modP\omega^n \equiv 1 \;mod\;P 的属于模 PP 的完全剩余系的非负整数 ω\omega 。nn 次单位根正好有 nn 个,它们是 g(P−1)knmodPg^{\frac{(P-1)k}{n}}\;mod\;P,k=0,1,⋯,n−1k=0,1,\cdots,n-1。
值 ωn≡gP−1nmodP\omega_n \equiv g^{\frac{P-1}{n}} \;mod\;P 称为主 nn 次单位根;所有的其他 nn 次单位根都是 ωn\omega_n 的幂。
同样,和复数域上的单位复根类似,nn 个 nn 次单位根ω0n,ω1n,⋯,ωn−1n\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1} 在模 PP 的乘法运算下形成了一个 nn 阶的循环群。
类比复数域上的单位复根,给出 nn 次单位根的重要性质。
相消引理:对任何整数 n≥0n \geq 0 , k≥0k \geq 0 , d>0d>0 , ωdkdn≡ωknmodP\omega_{dn}^{dk} \equiv \omega_{n}^{k} \;mod\;P ,由定义易证。
相消引理推论:ωn2n≡−1modP\omega_{n}^{\frac{n}{2}} \equiv -1\;mod\;P .
折半引理:如果 n>0n>0 为偶数,nn 个 nn 次单位根的平方等于 n2\frac{n}{2} 个 n2\frac{n}{2} 次单位根,根据相消引理有 (ωk+n2n)2≡ωk+n2n2≡ωkn2≡(ωkn)2modP(\omega_n^{k+\frac{n}{2}})^2\equiv\omega_{\frac{n}{2}}^{k+\frac{n}{2}}\equiv\omega_{\frac{n}{2}}^{k}\equiv(\omega_{n}^{k})^2\;mod\;P,其中 0≤k<n20\leq k 。
由相消引理的推论还能得出折半引理推论:ωkn≡−ωk+n2nmodP\omega_{n}^{k}\equiv-\omega_n^{k+\frac{n}{2}}\;mod\;P,其中 0≤k<n20\leq k 。
求和引理:对任意整数 n≥1n\geq1 和不能被 nn 整除的 kk ,有 ∑i=0n−1(ωkn)i≡0modP\sum\limits_{i=0}^{n-1}(\omega_{n}^{k})^i\equiv0\;mod\;P ,由等比数列求和公式易证。
求和引理推论:对任意整数 n≥1n\geq1 和整数 kk ,有 ∑i=0n−1(ωkn)i≡n[n|k]modP\sum\limits_{i=0}^{n-1}(\omega_{n}^{k})^i\equiv n[n|k]\;mod\;P ,其中 [n|k][n|k] 当且仅当 kk 能被 nn 整除时等于 11 ,否则等于 00。
NTT
回顾一下,我们希望计算次数界为 nn 的多项式 A(x)=∑i=0n−1aixiA(x)=\sum\limits_{i=0}^{n-1} a_i x^i 在 ω0n\omega_n^0,ω1n\omega_n^1,⋯\cdots,ωn−1n\omega_n^{n-1} 处的值。不是一般性地假定 nn 是 22 的幂,因为给定的次数界总可以增大;如果需要,总可以添加值为零的新的高阶系数。假定已知 AA 的系数形式:a=(a0,a1,⋯,an−1)a=(a_0,a_1,\cdots,a_{n-1}) 。对 k=0k=0,11,⋯\cdots,n−1n-1,定义结果 yky_k 如下:
yk≡A(ωkn)≡∑i=0n−1aiωiknmodPy_k\equiv A(\omega_n^k)\equiv \sum\limits_{i=0}^{n-1} a_i\omega_n^{ik} \;mod\;P
向量 y=(y0,y1,⋯,yn−1)y=(y_0,y_1,\cdots,y_{n-1}) 是系数向量 a=(a0,a1,⋯,an−1)a=(a_0,a_1,\cdots,a_{n-1}) 的数论变换( Number Theoretic Transform, NTT),也写作 y=NTTn(a)y=NTT_n(a)。
FNTT
通过一种称为快速数论变换(FNTT),就可以在 O(nlogn) O(nlogn) 的时间内计算出 NTTn(a)NTT_n(a),而直接的方法所需时间为 O(n2)O(n^2) 。FNTT主要是利用了单位根的性质。
FNTT方法运用了分治策略,它用 A(x)A(x) 中偶数下标的系数与奇数下标的系数,分别定义了两个新的次数界为 n2\frac{n}{2} 的多项式 A[0]A^{[0]} 和 A[1]A^{[1]} :
A[0](x)≡a0+a2x+a4x2+⋯+an−2xn2−1modPA^{[0]}(x)\equiv a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1} \;mod\;P
A[1](x)≡a1+a3x+a5x2+⋯+an−1xn2−1modPA^{[1]}(x)\equiv a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1} \;mod\;P
注意, A[0]A^{[0]} 包括 AA 中所有偶数下标的系数(下标的相应二进制数对应位置为 00), A[1]A^{[1]} 包括 AA 中所有奇数下标的系数(下标的相应二进制数对应位置为 11)。有下式:
A(x)≡A[0](x2)+xA[1](x2)modPA(x)\equiv A^{[0]}(x^2)+xA^{[1]}(x^2) \;mod\;P
这样,求 A(x)A(x) 在 ω0n\omega_n^0, ω1n\omega_n^1, ⋯\cdots, ωn−1n\omega_n^{n-1} 处的值的问题就转换为:
(1)(1) 求次数界为 n2\frac{n}{2} 的多项式 A[0]A^{[0]} 和 A[1]A^{[1]} 在 ω0n2\omega_{\frac{n}{2}}^0, ω1n2\omega_{\frac{n}{2}}^1, ⋯\cdots, ωn2−1n2\omega_{\frac{n}{2}}^{\frac{n}{2}-1} 处的值,然后
(2)(2) 根据下面的公式将 (1)(1) 的结果进行结合,对于 0≤k<n20\leq k 有:
A(ωkn)≡A[0](ωkn2)+ωknA[1](ωkn2)modPA(\omega_{n}^{k})\equiv A^{[0]}(\omega_{\frac{n}{2}}^{k})+\omega_n^kA^{[1]}(\omega_{\frac{n}{2}}^{k}) \;mod\;P
A(ωn2+kn)≡A[0](ωkn2)−ωknA[1](ωkn2)modPA(\omega_{n}^{\frac{n}{2}+k})\equiv A^{[0]}(\omega_{\frac{n}{2}}^{k})-\omega_n^kA^{[1]}(\omega_{\frac{n}{2}}^{k}) \;mod\;P
递归基为 n=1n=1,对于一个次数界等于 11 的多项式 AA ,有 A(x)=a0A(x)=a_0 ,其中不含有与 xx 有关的项,因此其在任意点处对应的点值都为 a0a_0 。
时间复杂度为 T(n)=2T(n2)+O(n)T(n)=2T(\frac{n}{2})+O(n),根据主定理有 T(n)=O(nlogn)T(n)=O(nlogn) (或者可以这样分析,递归中每一层的时间复杂度都是 O(n)O(n) ,一共有 lognlogn 层)。
因此,运用快速数论变换,可以在 O(nlogn)O(nlogn) 的时间内,求出次数界为 nn 的多项式在 nn 次单位根处的值。
对单位根进行插值
现在来说明如何在单位根处进行插值,以便把一个多项式从点值表示转化为系数表示。按如下方式进行插值:把NTT写成一个矩阵方程,然后再检查其逆矩阵的形式。
根据 yk≡A(ωkn)≡∑i=0n−1aiωiknmodPy_k\equiv A(\omega_n^k) \equiv \sum\limits_{i=0}^{n-1} a_i \omega_n^{ik} \;mod\;P ,其中 k=0k=0,11,⋯\cdots ,n−1n-1 ,可以把NTT写成矩阵积 y=Vnay=V_na,其中 VnV_n 是 ωn\omega_n 的适当幂组成的一个范德蒙德矩阵:
\begin{equation} \left[ \begin{matrix}y_0\\y_1\\y_2\\y_3\\\vdots\\y_{n-1} \end{matrix} \right] = \left[ \begin{matrix} 1&1&1&1&\cdots&1\\1&\omega_n&\omega_n^2&\omega_n^3&\cdots&\omega_n^{n-1}\\1&\omega_n^2&\omega_n^4&\omega_n^6&\cdots&\omega_n^{2(n-1)}\\1&\omega_n^3&\omega_n^6&\omega_n^9&\cdots&\omega_n^{3(n-1)}\\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\omega_n^{3(n-1)}&\cdots&\omega_n^{(n-1)(n-1)}\\ \end{matrix} \right] \left[ \begin{matrix} a_0\\a_1\\a_2\\a_3\\\vdots\\a_{n-1} \end{matrix} \right] \end{equation}
对 i,j=0,1,⋯,n−1i,j=0,1,\cdots,n-1, VnV_n 的 (i,j)(i,j) 处元素为 ωijn\omega_n^{ij},并且 VnV_n 中元素的幂形成一张乘法运算表。
对于逆运算 a=NTT−1n(y)a=NTT_n^{-1}(y),通过把 yy 乘以逆矩阵 V−1nV_n^{-1} 来进行处理。
定理:对 i,j=0,1,⋯,n−1i,j=0,1,\cdots,n-1,V−1nV_n^{-1} 的 (i,j)(i,j) 处的元素为 ω−ijnnmodP\frac{\omega_{n}^{-ij}}{n} \;mod\;P 。
证明: 对 i,j=0,1,⋯,n−1i,j=0,1,\cdots,n-1,V−1nVnV_n^{-1}V_n 的 (i,j)(i,j) 处的元素为 ∑k=0n−1ω−iknωkjnn=∑k=0n−1ωk(j−i)nn\sum\limits_{k=0}^{n-1} \frac{\omega_n^{-ik}\omega_n^{kj}}{n}=\sum\limits_{k=0}^{n-1}\frac{\omega_n^{k(j-i)}}{n} ,由求和引理推论有 ∑k=0n−1ωk(j−i)nn=[n|(j−i)]\sum\limits_{k=0}^{n-1}\frac{\omega_n^{k(j-i)}}{n}=[n|(j-i)],因此当且仅当 i=ji=j 时该式为1,否则为0,即 V−1nVn=EV_n^{-1}V_n=E,命题得证。
在求出逆矩阵 V−1nV_n^{-1} 之后, NTT−1n(y)NTT_n^{-1}(y) 可由下式给出:
ai≡1n∑j=0n−1yjω−ijnmodPa_i \equiv \frac{1}{n} \sum\limits_{j=0}^{n-1} y_j \omega_n^{-ij} \;mod\;P
其中 i=0,1,⋯,n−1i=0,1,\cdots,n-1 。
对比 NTT, INTT 只要把对应位置用到的单位根 ωn\omega_n 的幂次换成相反数(根据单位根的性质有 ωin≡ωimodnnmodP\omega_n^{i} \equiv \omega_n^{i\,mod\,n}\;mod\;P),最后再都乘上一个 nn 在模 PP 意义下的乘法逆元即可。
因此INTT的时间复杂度和NTT一样也是 O(nlogn)O(nlogn) 。
卷积定理:对任意两个长度为 nn 的向量 aa 和 bb,其中 nn 是 22 的幂,
a⨂b=NTT−12n(NTT2n(a)⋅NTT2n(b))a \bigotimes b = NTT^{-1}_{2n}(NTT_{2n}(a) \cdot NTT_{2n}(b))
其中 aa 和 bb 用 00 扩充使其长度达到 2n2n,⋅\cdot 表示 22 个 2n2n 个元素组成的向量的点乘。
因此,利用这种方法,可以在 O(nlogn)O(nlogn) 的时间内实现多项式乘法。
有效的FNTT实现
使用迭代的方式避免递归,预处理出原向量中每个位置在最底层被交换到的位置,预处理出 ωn\omega_n 的幂次方便进行蝴蝶操作。预处理完之后算法只需要先通过 O(n)O(n) 次交换得到最底层的情况,然后自底向上通过蝴蝶操作进行合并即可。
假设进行 FNTT 的多项式的次数界为 n=2kn=2^k ,则原本位置标号为 x=0,1,⋯,n−1x=0,1,\cdots,n-1 的元素在第 i=0,1,⋯,k−1i=0,1,\cdots,k-1 层会被分到 A[1]A^{[1]} 中的条件为 [x2i][\frac{x}{2^i}] 为奇数,因此对应的 xx 在最底层被交换到的位置 yx=∑i=0k−1([x2i]mod2)2k−1−iy_x=\sum\limits_{i=0}^{k-1} ([\frac{x}{2^i}]\;mod\;2)2^{k-1-i},因此 yxy_x 是 xx 在 kk 位二进制下反转后得到的数字。
通过递推求 yiy_i 的方法是 yi=[yi22]+[imod2]2k−1y_i=[\frac{y_{\frac{i}{2}}}{2}]+[i\;mod\;2]2^{k-1},i=1,2,⋯,n−1i=1,2,\cdots,n-1,y0=0y_0=0 。
探索类问题7(1)
给出两个多项式 ∑n=023n2xn\sum\limits_{n=0}^{2^3} n^2 x^n,∑n=024nxn\sum\limits_{n=0}^{2^4} n x^n,用快速傅里叶变换实现多项式相乘,并计算时间复杂度。
令 A(x)=∑i=023−1(i+1)2xiA(x)=\sum\limits_{i=0}^{2^3-1}(i+1)^2x^i,B(x)=∑i=024−1(i+1)xiB(x)=\sum\limits_{i=0}^{2^4-1} (i+1)x^i ,则原问题等价于求 x2A(x)B(x)x^2A(x)B(x),即计算次数界为 232^3 的多项式 AA 和次数界为 242^4 的多项式 BB 的乘积,然后将得到的系数向量依次向右移动 22 位,并令前两个元素为 00 即可。
令 m=25m=2^5,取 P=998244353=119×223+1P=998244353=119\times2^{23}+1,PP 为素数,最小的原根为 g=3g=3 ,则在模 PP 的意义下进行计算 a⨂b=NTT−1m(NTTm(a)⋅NTTm(b))a\bigotimes b=NTT^{-1}_{m}(NTT_m(a)\cdot NTT_m(b)),其中 aa 和 bb 分别是 AA 和 BB 扩充到 mm 的系数向量,由于 m4=220=1048576<Pm^4=2^{20}=1048576
,积多项式的每一项的系数都小于 PP ,取模没有影响,因此得到的系数就是原本的系数。
时间复杂度为 O(mlogm)。O(mlogm) 。
用来解决问题 (1)(1) 的代码如下:
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;#define max_N (1<<10)
#define P 998244353
#define g 3
//g is the Primitive Root of Pint n,inv_n;
//n is a power of 2 which represents the.degree-bound of the polynomial
//inv_n is the multiplicative inverse of n mod Pint brc[max_N];//binary reversing copy int w[2][max_N];//twiddle factors,w[0] for DFT while w[1] for IDFTint qpow(int x,int k){//an effective algorithm to calculate x^k mod Pint ans=1;for(;k;x=1ll*x*x%P,k>>=1)if(k&1)ans=1ll*ans*x%P;return ans;
}void fft_init(int l){//n is equal to 1<<l//calculating the multiplicative inverse of n mod Pinv_n=qpow(n,P-2);//due to Fermat Theory//calculate twiddle factors mod Pint G=qpow(g,(P-1)/n);w[0][0]=w[1][0]=1;for(int i=1;i<n;++i)w[0][i]=1ll*w[0][i-1]*G%P;for(int i=1;i<n;++i)w[1][i]=w[0][n-i];//calculate binary reversionsfor(int i=1;i<n;++i)brc[i]=(brc[i>>1]>>1)^((i&1)<<(l-1));}void fft(int*a,int t){//Fast Fourier Transform
//if t==0, do the Discrete Fourier Transform
//if t==1, do the Inverse Discrete Fourier Transform//binary reversing copyfor(int i=1;i<n;++i)if(i<brc[i])swap(a[i],a[brc[i]]);//the butterfly operationfor(int i=2;i<=n;i<<=1)for(int j=0;j<n;j+=i)for(int k=0;k<(i>>1);++k){int tmp=1ll*a[(i>>1)+j+k]*w[t][n/i*k]%P;a[(i>>1)+j+k]=(a[j+k]-tmp+P)%P;a[j+k]=(a[j+k]+tmp)%P;}//do not forget this step when doing IDFTif(t)for(int i=0;i<n;++i)a[i]=1ll*a[i]*inv_n%P;
}int a[max_N],b[max_N],c[max_N];int main(){freopen("answer1.txt","w",stdout);//output to a file named "answer1.txt"//initialize the coeffient of two polynomialsfor(int i=0;i<(1<<3);++i)a[i]=(i+1)*(i+1);for(int i=0;i<(1<<4);++i)b[i]=i+1;n=1<<5,fft_init(5);//initialize fftfft(a,0),fft(b,0);//compute the point-value representation of a and bfor(int i=0;i<n;++i)c[i]=1ll*a[i]*b[i]%P;//a*b mod P can be linearly calculated in point-value representationfft(c,1);//compute the coeffient of a*b//ascending print the coeffient//the answer is x^2*a*bputs("0\n0");for(int i=0;i<n;++i)printf("%d\n",c[i]);puts("");return 0;
}
该代码将积多项式的系数逐行升序输出到当前目录下的文件”Answer1.txt”中,文件共有 m+2m+2 行,第 i=1,2,⋯,m+1i=1,2,\cdots,m+1 行为积多项式中 xi−1x^{i-1} 的系数。
另附验证”Answer1.txt”的正确性的代码如下:
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;int gint(){char c; int f=1;while(c=getchar(),c<48||c>57)if(c=='-')f=-1;int x=0;for(;c>47&&c<58;c=getchar()){x=x*10+c-48;}return x*f;
}#define max_N 2333int n,a[max_N],b[max_N],c[max_N];int main(){freopen("answer1.txt","r",stdin);for(int i=0;i<(1<<3);++i)a[i]=(i+1)*(i+1);for(int i=0;i<(1<<4);++i)b[i]=i+1;n=1<<5;for(int i=0;i<n;++i)for(int j=0;j<=i;++j){c[i]+=a[j]*b[i-j];}gint(),gint();for(int i=0,x;i<n;++i)if(gint()!=c[i]){puts("WrongAnswer!");for(;;);}puts("Accepted!");for(;;);return 0;
}
该验证代码用 O(n2)O(n^2) 的卷积求出积函数的系数向量,并且从”Answer1.txt”中读入系数与该程序计算出的系数进行比较,如果发现不同,则输出”WrongAnswer!”,如果没有不同,则会输出”Accepted!”。
探索类问题7(2)
给定多项式 ∑n=0220n2xn\sum\limits_{n=0}^{2^{20}}n^2x^n,∑n=0220nxn\sum\limits_{n=0}^{2^{20}}nx^n,用快速傅里叶变换实现多项式相乘,并计算时间复杂度。
同 (1)(1) ,令 A(x)=∑i=0220−1(i+1)2xiA(x)=\sum\limits_{i=0}^{2^{20}-1}(i+1)^2x^i,B(x)=∑i=0220−1(i+1)xiB(x)=\sum\limits_{i=0}^{2^{20}-1}(i+1)x^i,则问题转换为求 x2A(x)B(x)x^2A(x)B(x) ,两个次数界都为 2202^{20} 的多项式 AA 和 BB 相乘。
令 m=221m=2^{21},与 (1)(1) 不同,由于系数可能达到 280≈1.2×10242^{80}\approx 1.2 \times 10^{24} ,若找一个大于 102410^{24} 数量级的素数 PP 来进行计算,中间结果可以很大,现有的整型数据类型难以支持,就算使用高精度也会很大程度上影响程序的运行效率,我们需要另辟蹊径。
中国剩余定理
若正整数 m1,m2,⋯,mkm_1,m_2,\cdots,m_k 两两互质,则同余方程组
x≡a1modm1x\equiv a_1\;mod\;m_1
x≡a2modm2x\equiv a_2\;mod\;m_2
⋯\cdots
x≡akmodmkx\equiv a_k\;mod\;m_k
有整数解,并且在模 M=∏i=1kmiM=\prod\limits_{i=1}^k m_i 意义下有唯一解, x≡∑i=1kai⋅Mmi⋅invi(Mmi)modMx\equiv\sum\limits_{i=1}^k a_i \cdot \frac{M}{m_i} \cdot inv_i(\frac{M}{m_i}) \;mod\;M,其中, invi(x)⋅x≡1modmiinv_i(x)\cdot x \equiv 1 \;mod\;m_i。
因此,我们可以分别计算出积多项式在若干个不同的素数模意义下的系数向量,只要保证这些模数的乘积大于 2802^{80},就可以使用中国剩余定理还原出本来的系数。
取 P0=469762049=7×226+1P_0=469762049=7 \times 2^{26} +1,P1=998244353=119×223+1P_1=998244353=119\times2^{23}+1,P2=1004535809=479×221+1P_2=1004535809=479\times2^{21}+1,这三个模数都是可以用于NTT的素数,并且最小的原根都为 g=3g=3 ,P2P_2 做加法刚好不会溢出 int类型,P=P0×P1×P2≈4.7×1026P=P_0\times P_1 \times P_2 \approx 4.7 \times 10^{26}。
分别计算出积多项式在模P0P_0,P1P_1,P2P_2 意义下的系数,再使用中国剩余定理还原真正的系数,由于 280<P<2127−1≈1.7×10382^{80}
,合并时可以使用64位系统自带的__int128数据类型进行计算,但是输出时由于系统不支持直接输出,需要转成十进制表示的字符串进行输出。
一共进行了 33 次多项式乘法,时间复杂度为 O(nlogn)O(nlogn),用中国剩余定理的部分时间复杂度为 O(n)O(n),算法总的时间复杂度为 O(nlogn)O(nlogn) 。
用来解决问题 (2)(2) 的代码如下:
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;#define max_N (1<<21)const int P[]={469762049,998244353,1004535809};
//these prime numbers have the same primitive root
#define g 3
//the primitive root of Pint n,l,inv_n;
//n is a power of 2 which represents the.degree-bound of the polynomial
//n is equal to 1<<l, inv_n is the multiplicative inverse of n mod Pint brc[max_N];//binary reversing copyint w[2][max_N];//twiddle factors, w[0] for DFT, w[1] for IDFTint qpow(int x,int k,int P){//an effective algorithm to computing x^k mod Pint ans=1;for(;k;x=1ll*x*x%P,k>>=1)if(k&1)ans=1ll*ans*x%P;return ans;
}void fft_init(int k){//calculating the multiplicative inverse of n mod P[k]inv_n=qpow(n,P[k]-2,P[k]);//due to Fermat Theory//computing twiddle factors mod P[k]int G=qpow(g,(P[k]-1)/n,P[k]);w[0][0]=w[1][0]=1;for(int i=1;i<n;++i)w[0][i]=1ll*w[0][i-1]*G%P[k];for(int i=1;i<n;++i)w[1][i]=w[0][n-i];
}void fft(int*a,int d,int t){//Fast Fourier Transform mod P[d]//if t==0, do the Discrete Fourier Transform //if t==1, do the Inverse Discrete Fourier Transformfor(int i=1;i<n;++i)if(i<brc[i])swap(a[i],a[brc[i]]);//butterfuly operationfor(int i=2;i<=n;i<<=1)for(int j=0;j<n;j+=i)for(int k=0;k<(i>>1);++k){int tmp=1ll*a[(i>>1)+j+k]*w[t][n/i*k]%P[d];a[(i>>1)+j+k]=(a[j+k]-tmp+P[d])%P[d];a[j+k]=(a[j+k]+tmp)%P[d];}//do not forget this step when doing IDFT;if(t)for(int i=0;i<n;++i)a[i]=1ll*a[i]*inv_n%P[d];
}int a[max_N],b[max_N],c[max_N];__int128 ans[max_N];inline void print(__int128 x){//system can not print __int128 directly ...if(!x){putchar('0');return;}static char c[40];int l=0;while(x)c[++l]=48+x%10,x/=10;while(l)putchar(c[l--]);putchar('\n');
}int main(){freopen("answer2.txt","w",stdout);//output to the file named "answer2.txt"n=1<<21,l=21;//computing binary reversionsfor(int i=1;i<n;++i)brc[i]=(brc[i>>1]>>1)^((i&1)<<(l-1));__int128 p=__int128(P[0])*P[1]*P[2];for(int k=0;k<3;++k){//initialize the coeffient mod P[k] of two polynomialsfor(int i=0;i<(1<<20);++i)a[i]=1ll*(i+1)*(i+1)%P[k];for(int i=(1<<20);i<(1<<21);++i)a[i]=0;for(int i=0;i<(1<<20);++i)b[i]=i+1;for(int i=(1<<20);i<(1<<21);++i)b[i]=0;fft_init(k);//initialize fftfft(a,k,0),fft(b,k,0);//compute the point-value representation of a and bfor(int i=0;i<n;++i)c[i]=1ll*a[i]*b[i]%P[k];//a*b mod P can be linearly calculated in point-value representationfft(c,k,1);//compute the coeffient of a*b mod P[k]//compute the coeffient of a*b with Chinese remainder theorem__int128 N=1ll*P[(k+1)%3]*P[(k+2)%3];__int128 M=qpow(N%P[k],P[k]-2,P[k]);for(int i=0;i<n;++i)ans[i]=ans[i]+N*M*c[i]%p;}//ascending print the coeffient//the answer is x^2*a*bputs("0\n0");for(int i=0;i<n;++i)print(ans[i]%p);puts("");return 0;
}
该代码将积多项式的系数逐行升序输出到当前目录下的文件”Answer2.txt”中,文件共有 m+2m+2 行,第 i=1,2,⋯,m+1i=1,2,\cdots,m+1 行为积多项式中 xi−1x^{i-1} 的系数。
另附验证”Answer2.txt”的正确性的代码如下:
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;#define max_N (1<<21)+5
const int P[]={469762049,998244353,1004535809};int qpow(int x,int k,int P){//an effective algorithm to computing x^k mod Pint ans=1;for(;k;x=1ll*x*x%P,k>>=1)if(k&1)ans=1ll*ans*x%P;return ans;
}int n,c[3][max_N];int sum2(int n,int inv,int P){//sum i^2 mod P, inv=1/6 mod Preturn 1ll*n*(n+1)%P*(n<<1|1)%P*inv%P;
}int sum3(int n,int P){//sum i^3 mod Pint tmp=1ll*n*(n+1)/2%P;return 1ll*tmp*tmp%P;
}int main(){freopen("answer2.txt","r",stdin);n=1<<21;for(int i=0;i<n+2;++i){char ch;while(ch=getchar(),ch<48||ch>57);for(;ch>47&&ch<58;ch=getchar())for(int k=0;k<3;++k){c[k][i]=(1ll*c[k][i]*10+ch-48)%P[k];}}for(int k=0;k<3;++k){int inv=qpow(12,P[k]-2,P[k]);// 1/12int inv1=qpow(6,P[k]-2,P[k]);// 1/6int tmp1=sum2(1<<20,inv1,P[k]);int tmp2=sum3(1<<20,P[k]);for(int i=2,ans;i<n+2;++i){if(i<=(1<<20)){ans=1ll*i*i%P[k];ans=1ll*ans*(ans-1)%P[k];ans=1ll*ans*inv%P[k];ans=(ans+P[k])%P[k];}else{ans=1ll*i*(tmp1-sum2(i-(1<<20)-1,inv1,P[k]))%P[k];ans=(ans-tmp2+sum3(i-(1<<20)-1,P[k]))%P[k];ans=(ans+P[k])%P[k];}if(ans!=c[k][i]){printf("%d %d\n",k,i);printf("%d %d\n",ans,c[k][i]);puts("WrongAnswer!");for(;;);}}}puts("Accepted!");for(;;);return 0;
}
通过推公式的方法,用 O(n)O(n) 的时间计算出积多项式在P0P_0,P1P_1,P1P_1 模意义下的系数,从”Answer2.txt”中读入系数并对P0P_0,P1P_1,P1P_1 分别取模与该程序计算出的对应模意义下的系数进行比较,如果发现不同,则输出”WrongAnswer!”,如果没有不同,则会输出”Accepted!”,由中国剩余定理,这种验证方法的可靠性显然。
上面所提到的公式为:
ck=k2(k2−1)12c_k = \frac{k^2(k^2-1)}{12} ,k=0,1,⋯,220k=0,1,\cdots,2^{20}
ck=k(s2(220)−s2(k−220−1))−(s3(220)−s3(k−220−1))c_k = k(s_2(2^{20})-s_2(k-2^{20}-1)) - (s_3(2^{20})-s_3(k-2^{20}-1)),k=220+1,220+2,⋯,221−1k=2^{20}+1,2^{20}+2,\cdots,2^{21}-1
其中 s2(n)=∑i=0ni2=n(n+1)(2n+1)6s_2(n)=\sum\limits_{i=0}^{n} i^2=\frac{n(n+1)(2n+1)}{6}, s3(n)=∑i=0ni3=[n(n+1)2]2s_3(n)=\sum\limits_{i=0}^{n} i^3=\left[\frac{n(n+1)}{2}\right]^2 。
工科数学分析大作业(三) 傅里叶级数相关推荐
- web前端设计与开发大作业(三)----注册与登陆界面
一.登陆界面 <!DOCTYPE html> <html lang="en"> <head><meta charset="UTF ...
- 清华贵系的期末大作业:奋战三周,造台计算机!
大数据文摘授权转载自AI科技评论 作者 | 蒋宝尚 编辑丨陈彩娴 本科大三,正在学习计算机组成原理,能做个什么项目? 清华大学贵系说:造台计算机吧! 清华有门本科三年级必修课,名为<计算机组成原 ...
- 201609计算机控制技术作业三,计算机控制技术大作业2015..doc
计算机控制技术大作业2015. 深圳大学考试答题纸 (以论文.报告等形式考核专用)二○一四 -二○一五 学年度第 2 学期 课程编号1700470001课程名称计算机控制技术主讲教师评分学 号姓名专业 ...
- 第三次大作业-作业准备
在第二次大作业结束之后,第三次大作业-白盒测试实践在三周内进行. 由于大家刚进行完一次大作业的辛苦工作,于是本组决定第一周交由大家复习学习之前的黑盒测试实践.搜集资料.自发配置工作环境和工具.进行自己 ...
- C++面向对象程序设计大作业:魔兽世界(三):开战
C++面向对象程序设计大作业:魔兽世界(三):开战 问题描述 问题分析 代码 问题描述 问题来自于北京大学郭炜老师的C++慕课的大作业 魔兽世界的西面是红魔军的司令部,东面是蓝魔军的司令部.两个司令部 ...
- 密码学大作业(共三次)
第一次大作业 Many Time Pad #密文为十六进制字符串,应该先将其处理 ciphertexts = ["315c4eeaa8b5f8aaf9174145bf43e1784b8fa0 ...
- 高级软件工程第三次大作业(周帅)
一如软工深似海,从此硬件是故人 课前疑问的自答 课前,我的主要疑问是通过这门课对我的课题有什么帮助?我了解软件开发流程对于我的课题有什么意义?伴随着这门课的学习,我有了一个契机去主动了解Pycharm ...
- 《spark技术应用》课程期末考试大作业报告,使用eclipse完成求top值、文件排序、二次排序三个程序的个性化开发。
目录 一.选题的目的及要求... 4 二.设计思路... 4 三.主要内容及关键技术.. 5 四.制作步骤... 5 1.准备工作... 5 1.1在VMware中安装一台Ubuntu64位系 ...
- 魔兽世界(三):开战 北京大学慕课C++大作业 代码与思路
类的组织 总体思路: 自底向上为:Weapon -> Warrior -> HeadQuarter | City -> main函数 每个司令部含有一个Warrior的vector, ...
- HTML5期末大作业:酒店主题——绿色温泉度假酒店网页设计(8页) HTML+CSS+JavaScript 大三网页设计源码
HTML5期末大作业:酒店主题--绿色温泉度假酒店网页设计(8页) HTML+CSS+JavaScript 期末作业HTML代码 学生网页课程设计期末作业下载 web网页设计制作成品 常见网页设计作业 ...
最新文章
- Java中stringBuilder
- pelee yuface 手势模型
- Spring第一讲:初步了解Spring
- 如何解析lvx文档_建站零基础入门:手把手教你如何自助建站
- 干货:产品经理怎么做才能在需求评审中少挨打?
- 计算机应用技术环境评估,计算机应用教程(第7版)(Windows 7与Office 2007环境)习题解答与上机练习...
- 分析log及校准学习总结
- 车借给朋友好几次,满油的车每次还回来都是没油了,我觉得心里有些不舒服是我太计较吗?
- SQL Server 查询处理中的各个阶段
- android 项目 功能 源码 eclipse的
- zookeeper入门学习《一》
- POJ 3597 Polygon Division (DP)
- win10cmd重置系统_win10命令提示符一键还原修复系统
- Steam 游戏服务器无法连接 steam 游戏无法启动 打开 microsoft store 错误代码 0x80131500
- Springboot JUnit5 Controller 单元测试
- Lottie 动画在项目中的使用总结
- 使用Python绘制二元函数图像
- 计算机黑屏时间,电脑开机黑屏时间长怎么办?Win10开机黑屏时间很久的解决方法...
- 【数据集划分】误用shuffle,导致训练集和测试集掺混
- 阿里云的NoSQL存储服务OTS的应用分析
热门文章
- 食品进销存十大品牌排行榜新鲜出炉,来看看哪个最适合你
- 2022抖音日活用户超8亿,旅游商家如何从抖音获客?
- 思维导图模板创意可爱简单,模板资源分享
- 阿里异构离线数据同步工具/平台DataX
- 黑金花大理石_不同产地的黑金花大理石有哪些特点?
- 无监督学习——非负矩阵分解(NMF)
- 计算机二级不看教材只刷题可以吗,中级会计可以只看轻松过关不看教材吗
- 对两个等长升序的序列查找中位数
- Canvas API - 江苏黑马 - 博客园
- 服务器图标怎么显示在任务栏,win10任务栏右侧的图标如何显示或隐藏起来?附图文教程...