工科数学分析大作业(三) 傅里叶级数

研究课题：探索类问题7

给定多项式：(1)∑n=023n2xn,∑n=024nxn ; (2)∑n=0220n2xn,∑n=0220nxn ;
利用快速傅里叶变换实现多项式相乘，并分析算法的复杂度。

算法介绍：快速数论变换

模 P 意义下的多项式

A(x)=∑i=0n−1aixi，称 a0 , a1 , ⋯ , an−1 为多项式的系数。所有系数都属于模 P 的完全剩余系。如果多项式 A(x) 的最高次的非零系数为 ak ，则称 A(x) 的次数是 k . 任何一个严格大于一个多项式次数的整数都是这个多项式的次数界。

多项式乘法

如果多项式 A(x) 和 B(x) 都是次数界为 n 的多项式，则说它们的乘积是一个次数界为 2n−1 的多项式积 C(x)，并满足对所有属于定义域的 x，都有 C(x)=A(x)⋅B(x)modP .

多项式的表示方法

系数表示法：
对一个次数界为 n 的多项式 A(x)=∑i=0n−1aixi 来说，其系数表示法就是一个由系数组成的向量 a=(a0,a1,⋯,an−1) 。
系数表示法下的多项式乘法
假设 A 和 B 都是次数界为 n 的多项式，a 和 b 分别是 A 和 B 的系数向量，则 A 和 B 的乘积 C 的系数向量 c=a⨂b=(c0,c1,⋯.cn−1)，其中 ci≡∑j=0iaj⋅bi−jmodP 。根据这个定义计算两个次数界为 n 的多项式乘法的时间复杂度为 O(n2) 。
点值表示法：
一个次数界为 n 的多项式 A(x) 的点值表示就是 n 个点值对所形成的集合：{(x0,y0),(x1,y1),⋯,(xn−1,yn−1)}，其中所有 xk 互不相同，并且当 k=0,1,⋯,n−1 时有yk=A(xk)modP。一个多项式可以有很多不同的点值表示。
如果任意选取 {xk}，根据霍纳法则（秦九韶算法）将 A(x) 从系数表示法转换为点值表示法（求值）的时间复杂度为 O(n2) 。如果巧妙的选取 {xk}，可以将转换的时间复杂度降为 O(nlogn) 。
求值计算的逆（从一个多项式的点值表示确定其系数表示中的系数）称为插值。

多项式插值的唯一性

对于任意 n 个相异点值对组成的集合 {(x0,y0),(x1,y1),⋯,(xn−1,yn−1)}，存在唯一的次数界为 n 的多项式 A(x)，满足 yk≡A(xk)modP，k=0,1,⋯,n−1 。
证明： 证明过程的基础是某个矩阵存在逆矩阵，yk=A(xk)modP，k=0,1,⋯,n−1 等价于矩阵方程

⎡⎣⎢⎢⎢⎢⎢11⋮1x0x1⋮xn−1x20x21⋮x2n−1⋯⋯⋱⋯xn−10xn−11⋮xn−1n−1⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢y0y1⋮yn−1⎤⎦⎥⎥⎥⎥⎥
左边的矩阵是范德蒙德矩阵，表示为V(x0,x1,⋯,xn−1)。该矩阵的行列式的值为 ∏0≤i<j≤n−1(xj−xi)，因此，如果 xk 相异，则该矩阵是可逆的。因此，对给定的唯一的点值表示，就能够求出系数 ai 的解：a=V(x0,x1,⋯,xn−1)−1y 。

点值表示法下的多项式乘法

点值表示法对于多项式乘法是很方便的。如果 C(x)=A(x)B(x)，则对任意点 xk，有 C(xk)=A(xk)B(xk)，并且对 A 的点值表示和 B 的点值表示进行逐点相乘，就可以得到 C 的点值表示。假设 A 和 B 的次数界都为 n，则 C 的次数界为 2n，因此要得到 C 的点值表示需要 2n 个点值对。因此，必须对 A 和 B 的点值表示进行“扩充”，使每个多项式都包含 2n 个点值对。
如果已知 A 的扩充点值表示：{(x0,y0),(x1,y1),⋯,(x2n−1,y2n−1)} 和 B 的相应扩充点值表示：{(x0,y′0),(x1,y′1),⋯,(x2n−1,y′2n−1)}，则 C 的点值表示为：{(x0,y0y′0modP),(x1,y1y′1modP),⋯,(x2n−1,y2n−1y′2n−1modP)} 。
如果已知两个扩充点值形式的输入多项式，使其相乘而得到点值形式的结果需要 O(n) 的时间。

原根

设 m 是正整数，a 是整数，若 a 模 m 的阶等于 φ(m) ，则称 a 为模 m 的一个原根（其中 φ(m) 是 m 的欧拉函数）。
假设 g 是 P 的一个原根，则 g0,g,g2,⋯,gP−2在模 P 的乘法运算下形成了一个 P−1 阶循环群，其中 P 为素数，且有 gP−1≡g0≡1modP 。

系数形式表示的多项式的快速乘法

能否采用利用关于点值形式的表示的多项式的线性时间乘法算法，来加快系数形式表示的多项式乘法运算的速度呢？答案依赖于能否快速把一个多项式在系数形式和点值形式之间转换。
当模数特殊的时候，当 P=r2k+1，且 P 为素数，存在原根 g 的时候，如果选取 gP−1n 作为单位根，则可以把一个次数界为 n 的多项式在两种表示法之间转化所需的时间压缩为 O(nlogn) （这里要求 n=2j 并且 j≤k）。
本文中之后所提到的所有模数 P 均默认为满足上述条件。

NTT和FNTT

模P意义下的单位根

类似于复数域上单位复根的定义，在模 P 的意义下定义 n 次单位根为满足 ωn≡1modP 的属于模 P 的完全剩余系的非负整数 ω 。n 次单位根正好有 n 个，它们是 g(P−1)knmodP，k=0,1,⋯,n−1。
值 ωn≡gP−1nmodP 称为主 n 次单位根；所有的其他 n 次单位根都是 ωn 的幂。
同样，和复数域上的单位复根类似，n 个 n 次单位根ω0n,ω1n,⋯,ωn−1n 在模 P 的乘法运算下形成了一个 n 阶的循环群。

类比复数域上的单位复根，给出 n 次单位根的重要性质。

相消引理：对任何整数 n≥0 ， k≥0 ， d>0 ， ωdkdn≡ωknmodP ，由定义易证。

相消引理推论：ωn2n≡−1modP .

折半引理：如果 n>0 为偶数，n 个 n 次单位根的平方等于 n2 个 n2 次单位根，根据相消引理有 (ωk+n2n)2≡ωk+n2n2≡ωkn2≡(ωkn)2modP，其中 0≤k<n2 。

由相消引理的推论还能得出折半引理推论：ωkn≡−ωk+n2nmodP，其中 0≤k<n2 。

求和引理：对任意整数 n≥1 和不能被 n 整除的 k ，有 ∑i=0n−1(ωkn)i≡0modP ，由等比数列求和公式易证。

求和引理推论：对任意整数 n≥1 和整数 k ，有 ∑i=0n−1(ωkn)i≡n[n|k]modP ，其中 [n|k] 当且仅当 k 能被 n 整除时等于 1 ，否则等于 0。

NTT
回顾一下，我们希望计算次数界为 n 的多项式 A(x)=∑i=0n−1aixi 在 ω0n，ω1n，⋯，ωn−1n 处的值。不是一般性地假定 n 是 2 的幂，因为给定的次数界总可以增大；如果需要，总可以添加值为零的新的高阶系数。假定已知 A 的系数形式：a=(a0,a1,⋯,an−1) 。对 k=0，1，⋯，n−1，定义结果 yk 如下：

yk≡A(ωkn)≡∑i=0n−1aiωiknmodP
向量 y=(y0,y1,⋯,yn−1) 是系数向量 a=(a0,a1,⋯,an−1) 的数论变换（Number Theoretic Transform, NTT），也写作 y=NTTn(a)。

FNTT
通过一种称为快速数论变换（FNTT），就可以在 O(nlogn) 的时间内计算出 NTTn(a)，而直接的方法所需时间为 O(n2) 。FNTT主要是利用了单位根的性质。
FNTT方法运用了分治策略，它用 A(x) 中偶数下标的系数与奇数下标的系数，分别定义了两个新的次数界为 n2 的多项式 A[0] 和 A[1] ：

A[0](x)≡a0+a2x+a4x2+⋯+an−2xn2−1modP
A[1](x)≡a1+a3x+a5x2+⋯+an−1xn2−1modP
注意，A[0] 包括 A 中所有偶数下标的系数（下标的相应二进制数对应位置为 0），A[1] 包括 A 中所有奇数下标的系数（下标的相应二进制数对应位置为 1）。有下式：
A(x)≡A[0](x2)+xA[1](x2)modP
这样，求 A(x) 在 ω0n，ω1n，⋯，ωn−1n 处的值的问题就转换为：
(1) 求次数界为 n2 的多项式 A[0] 和 A[1] 在 ω0n2，ω1n2，⋯，ωn2−1n2 处的值，然后
(2) 根据下面的公式将 (1) 的结果进行结合，对于 0≤k<n2 有：
A(ωkn)≡A[0](ωkn2)+ωknA[1](ωkn2)modP
A(ωn2+kn)≡A[0](ωkn2)−ωknA[1](ωkn2)modP
递归基为 n=1，对于一个次数界等于 1 的多项式 A ，有 A(x)=a0 ，其中不含有与 x 有关的项，因此其在任意点处对应的点值都为 a0 。

时间复杂度为 T(n)=2T(n2)+O(n)，根据主定理有 T(n)=O(nlogn) (或者可以这样分析，递归中每一层的时间复杂度都是 O(n) ，一共有 logn 层)。

因此，运用快速数论变换，可以在 O(nlogn) 的时间内，求出次数界为 n 的多项式在 n 次单位根处的值。

对单位根进行插值
现在来说明如何在单位根处进行插值，以便把一个多项式从点值表示转化为系数表示。按如下方式进行插值：把NTT写成一个矩阵方程，然后再检查其逆矩阵的形式。
根据 yk≡A(ωkn)≡∑i=0n−1aiωiknmodP ，其中 k=0，1，⋯ ，n−1 ，可以把NTT写成矩阵积 y=Vna，其中 Vn 是 ωn 的适当幂组成的一个范德蒙德矩阵：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢y0y1y2y3⋮yn−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1111⋮11ωnω2nω3n⋮ωn−1n1ω2nω4nω6n⋮ω2(n−1)n1ω3nω6nω9n⋮ω3(n−1)n⋯⋯⋯⋯⋱⋯1ωn−1nω2(n−1)nω3(n−1)n⋮ω(n−1)(n−1)n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢a0a1a2a3⋮an−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥

对 i,j=0,1,⋯,n−1，Vn 的 (i,j) 处元素为 ωijn，并且 Vn 中元素的幂形成一张乘法运算表。
对于逆运算 a=NTT−1n(y)，通过把 y 乘以逆矩阵 V−1n 来进行处理。

定理：对 i,j=0,1,⋯,n−1，V−1n 的 (i,j) 处的元素为 ω−ijnnmodP 。
证明： 对 i,j=0,1,⋯,n−1，V−1nVn 的 (i,j) 处的元素为 ∑k=0n−1ω−iknωkjnn=∑k=0n−1ωk(j−i)nn，由求和引理推论有 ∑k=0n−1ωk(j−i)nn=[n|(j−i)]，因此当且仅当 i=j 时该式为1，否则为0，即 V−1nVn=E，命题得证。

在求出逆矩阵 V−1n 之后， NTT−1n(y) 可由下式给出：

ai≡1n∑j=0n−1yjω−ijnmodP
其中 i=0,1,⋯,n−1 。

对比 NTT， INTT 只要把对应位置用到的单位根 ωn 的幂次换成相反数（根据单位根的性质有 ωin≡ωimodnnmodP），最后再都乘上一个 n 在模 P 意义下的乘法逆元即可。
因此INTT的时间复杂度和NTT一样也是 O(nlogn) 。

卷积定理：对任意两个长度为 n 的向量 a 和 b，其中 n 是 2 的幂，

a⨂b=NTT−12n(NTT2n(a)⋅NTT2n(b))
其中 a 和 b 用 0 扩充使其长度达到 2n，⋅ 表示 2 个 2n 个元素组成的向量的点乘。

因此，利用这种方法，可以在 O(nlogn) 的时间内实现多项式乘法。

有效的FNTT实现

使用迭代的方式避免递归，预处理出原向量中每个位置在最底层被交换到的位置，预处理出 ωn 的幂次方便进行蝴蝶操作。预处理完之后算法只需要先通过 O(n) 次交换得到最底层的情况，然后自底向上通过蝴蝶操作进行合并即可。

假设进行 FNTT 的多项式的次数界为 n=2k ，则原本位置标号为 x=0,1,⋯,n−1 的元素在第 i=0,1,⋯,k−1 层会被分到 A[1] 中的条件为 [x2i] 为奇数，因此对应的 x 在最底层被交换到的位置 yx=∑i=0k−1([x2i]mod2)2k−1−i，因此 yx 是 x 在 k 位二进制下反转后得到的数字。

通过递推求 yi 的方法是 yi=[yi22]+[imod2]2k−1，i=1,2,⋯,n−1，y0=0 。

探索类问题7(1)

给出两个多项式 ∑n=023n2xn，∑n=024nxn，用快速傅里叶变换实现多项式相乘，并计算时间复杂度。

令 A(x)=∑i=023−1(i+1)2xi，B(x)=∑i=024−1(i+1)xi ，则原问题等价于求 x2A(x)B(x)，即计算次数界为 23 的多项式 A 和次数界为 24 的多项式 B 的乘积，然后将得到的系数向量依次向右移动 2 位，并令前两个元素为 0 即可。

令 m=25，取 P=998244353=119×223+1，P 为素数，最小的原根为 g=3 ，则在模 P 的意义下进行计算 a⨂b=NTT−1m(NTTm(a)⋅NTTm(b))，其中 a 和 b 分别是 A 和 B 扩充到 m 的系数向量，由于 m4=220=1048576<P，积多项式的每一项的系数都小于 P ，取模没有影响，因此得到的系数就是原本的系数。

时间复杂度为 O(mlogm)。

用来解决问题 (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 P    int 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 P    inv_n=qpow(n,P-2);//due to Fermat Theory    //calculate twiddle factors mod P    int 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 reversions    for(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 copy    for(int i=1;i<n;++i)if(i<brc[i])swap(a[i],a[brc[i]]);    //the butterfly operation    for(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 IDFT    if(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 polynomials    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,fft_init(5);//initialize fft    fft(a,0),fft(b,0);    //compute the point-value representation of a and b    for(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 representation    fft(c,1);    //compute the coeffient of a*b    //ascending print the coeffient    //the answer is x^2*a*b    puts("0\n0");    for(int i=0;i<n;++i)printf("%d\n",c[i]);    puts("");    return 0;}

该代码将积多项式的系数逐行升序输出到当前目录下的文件”Answer1.txt”中，文件共有 m+2 行，第 i=1,2,⋯,m+1 行为积多项式中 xi−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) 的卷积求出积函数的系数向量，并且从”Answer1.txt”中读入系数与该程序计算出的系数进行比较，如果发现不同，则输出”WrongAnswer!”，如果没有不同，则会输出”Accepted!”。

探索类问题7(2)

给定多项式 ∑n=0220n2xn，∑n=0220nxn，用快速傅里叶变换实现多项式相乘，并计算时间复杂度。

同 (1) ，令 A(x)=∑i=0220−1(i+1)2xi，B(x)=∑i=0220−1(i+1)xi，则问题转换为求 x2A(x)B(x) ，两个次数界都为 220 的多项式 A 和 B 相乘。

令 m=221，与 (1) 不同，由于系数可能达到 280≈1.2×1024 ，若找一个大于 1024 数量级的素数 P 来进行计算，中间结果可以很大，现有的整型数据类型难以支持，就算使用高精度也会很大程度上影响程序的运行效率，我们需要另辟蹊径。

中国剩余定理

若正整数 m1,m2,⋯,mk 两两互质，则同余方程组

x≡a1modm1
x≡a2modm2
⋯
x≡akmodmk
有整数解，并且在模 M=∏i=1kmi 意义下有唯一解，x≡∑i=1kai⋅Mmi⋅invi(Mmi)modM，其中，invi(x)⋅x≡1modmi。

因此，我们可以分别计算出积多项式在若干个不同的素数模意义下的系数向量，只要保证这些模数的乘积大于 280，就可以使用中国剩余定理还原出本来的系数。

取 P0=469762049=7×226+1，P1=998244353=119×223+1，P2=1004535809=479×221+1，这三个模数都是可以用于NTT的素数，并且最小的原根都为 g=3 ，P2 做加法刚好不会溢出 int类型,P=P0×P1×P2≈4.7×1026。

分别计算出积多项式在模P0，P1，P2 意义下的系数，再使用中国剩余定理还原真正的系数，由于 280<P<2127−1≈1.7×1038，合并时可以使用64位系统自带的__int128数据类型进行计算，但是输出时由于系统不支持直接输出，需要转成十进制表示的字符串进行输出。

一共进行了 3 次多项式乘法，时间复杂度为 O(nlogn)，用中国剩余定理的部分时间复杂度为 O(n)，算法总的时间复杂度为 O(nlogn) 。

用来解决问题 (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 P    int 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 Transform    for(int i=1;i<n;++i)if(i<brc[i])swap(a[i],a[brc[i]]);    //butterfuly operation    for(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 reversions    for(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 polynomials        for(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 fft        fft(a,k,0),fft(b,k,0);        //compute the point-value representation of a and b        for(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 representation        fft(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*b    puts("0\n0");    for(int i=0;i<n;++i)print(ans[i]%p);    puts("");    return 0;}

该代码将积多项式的系数逐行升序输出到当前目录下的文件”Answer2.txt”中，文件共有 m+2 行，第 i=1,2,⋯,m+1 行为积多项式中 xi−1 的系数。

另附验证”Answer2.txt”的正确性的代码如下：

#include <stdio.h>#include <string.h>#include <algorithm>using namespace std;#define max_N (1<<21)+5const int P[]={469762049,998244353,1004535809};int qpow(int x,int k,int P){//an effective algorithm to computing x^k mod P    int 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 P    return 1ll*n*(n+1)%P*(n<<1|1)%P*inv%P;}int sum3(int n,int P){//sum i^3 mod P    int 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/12        int inv1=qpow(6,P[k]-2,P[k]);// 1/6        int 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) 的时间计算出积多项式在P0，P1，P1 模意义下的系数，从”Answer2.txt”中读入系数并对P0，P1，P1 分别取模与该程序计算出的对应模意义下的系数进行比较，如果发现不同，则输出”WrongAnswer!”，如果没有不同，则会输出”Accepted!”，由中国剩余定理，这种验证方法的可靠性显然。

上面所提到的公式为：

ck=k2(k2−1)12，k=0,1,⋯,220
ck=k(s2(220)−s2(k−220−1))−(s3(220)−s3(k−220−1))，k=220+1,220+2,⋯,221−1
其中 s2(n)=∑i=0ni2=n(n+1)(2n+1)6，s3(n)=∑i=0ni3=[n(n+1)2]2 。