多项式回归
 
多项式回归，回归函数是回归变量多项式的回归。多项式回归模型是线性回归模型的一种，此时回归函数关于回归系数是线性的。由于任一函数都可以用多项式逼近，因此多项式回归有着广泛应用。
 
直线回归研究的是一个因变量与一个自变量之间的回归问题，但在实际情况中，影响因变量的自变量往往不止一个，例如：羊毛的产量受到绵羊体重、体长、胸围等影响，因此需要进行一个因变量与多个自变量间的回归分析，即多元回归分析。
 
研究一个因变量与一个或多个自变量间多项式的回归分析方法，称为多项式回归（Polynomial Regression）。如果自变量只有一个时，称为一元多项式回归；如果自变量有多个时，称为多元多项式回归。在一元回归分析中，如果依变量y与自变量x的关系为非线性的，但是又找不到适当的函数曲线来拟合，则可以采用一元多项式回归。
 
一元m次多项式回归方程为：
 
二元二次多项式回归方程为：
 
 
 
 
上图的数据，我们可以使用一元2次多项式来拟合，首先，一个标准的一元高阶多项式函数如下：
 
 
m 表示多项式的阶数，表示 x 的 j 次幂，w 则代表该多项式的系数。
 
当我们使用上面的多项式去拟合散点时，需要确定两个要素，分别是：多项式系数 w 以及多项式阶数 m，这也是多项式的两个基本要素。当然也可以手动指定多项式的阶数m的大小，这样就只需要确定系数w的值了。得到以下公式：
 
 
如何求解该公式呢？这里使用Scipy的方法。
 
使用Scipy提供的最小二乘法函数得到最佳拟合参数：
 
该方法：最小化一组方程的平方和（即可以用来实现最小二乘法） 
 
import numpy as np
from scipy.optimize import leastsq

# 拟合数据集
x = [4, 8, 12, 25, 32, 43, 58, 63, 69, 79]
y = [20, 33, 50, 56, 42, 31, 33, 46, 65, 75]

def fun(p, x):
    """
    定义想要拟合的函数
    """
    w0,w1 = p  #从参数p获得拟合的参数
    # 如果是二次多项式则：w0,w1,w2 = p ；return w0 + w1*x + w2*x*x  以此类推
    return w0 + w1*x

def err(p, x, y):
    """
    定义误差函数
    """
    return fun(p,x) -y

#定义起始的参数 即从 y = 1*x+1 开始，其实这个值可以随便设，只不过会影响到找到最优解的时间
p0 = [1,1]   #p0 = [1,1,1]    w系数的个数[w0,w1,w2...]

#将list类型转换为 numpy.ndarray 类型，最初我直接使用
#list 类型,结果 leastsq函数报错，后来在别的blog上看到了，原来要将类型转
#换为numpy的类型

x1 = np.array(x)  
y1 = np.array(y)

xishu = leastsq(err, p0, args=(x1,y1))
print(xishu[0])
# xishu[0]，即为获得的参数
 
一般只要指定前三个参数就可以：
 
 
func 是我们自己定义的一个计算误差的函数，
 
 
x0 是计算的初始参数值
 
 
args 是指定func的其他参数
 
通过实践后观察，上面实现1次多项式拟合（2次多项式，p0则需要3个值）但其效果都不是很好，所以下面修改代码尝试N(大于2)项拟合
 
"""
	实现N次多项式拟合
"""
def fit_func(p, x):
	"""根据公式，定义 n 次多项式函数
	"""
	f = np.poly1d(p)   # 这里的np.poly1d函数是用来构造多项式使用的，默认格式为:ax**2+bx+c等，如：ax**3 + bx**2 + cx + d 以此类推
	return f(x)

def err_func(p, x, y):
	"""残差函数（观测值与拟合值之间的差距）
	"""
	ret = fit_func(p, x) - y
	return ret

def n_poly(n):
	"""n 次多项式拟合
	"""
	p_init = np.random.randn(n) # 生成 n 个随机数
	parameters = leastsq(err_func, p_init, args=(np.array(x), np.array(y)))
	return parameters[0]	# 返回多项式系数w0、w1、w2...

k = n_poly(3)	# 与上面的二次多项式结果一致，只是公式顺序不同

"""绘制出 3，4，5，6，7, 8 次多项式的拟合图像
"""
# 绘制拟合图像时需要的临时点
x_temp = np.linspace(0, 80, 10000)
# 绘制子图
fig, axes = plt.subplots(2, 3, figsize=(15,10))

axes[0,0].plot(x_temp, fit_func(n_poly(4), x_temp), 'r')
axes[0,0].scatter(x, y)
axes[0,0].set_title("m = 3")

axes[0,1].plot(x_temp, fit_func(n_poly(5), x_temp), 'r')
axes[0,1].scatter(x, y)
axes[0,1].set_title("m = 4")

axes[0,2].plot(x_temp, fit_func(n_poly(6), x_temp), 'r')
axes[0,2].scatter(x, y)
axes[0,2].set_title("m = 5")

axes[1,0].plot(x_temp, fit_func(n_poly(7), x_temp), 'r')
axes[1,0].scatter(x, y)
axes[1,0].set_title("m = 6")

axes[1,1].plot(x_temp, fit_func(n_poly(8), x_temp), 'r')
axes[1,1].scatter(x, y)
axes[1,1].set_title("m = 7")

axes[1,2].plot(x_temp, fit_func(n_poly(9), x_temp), 'r')
axes[1,2].scatter(x, y)
axes[1,2].set_title("m = 8")
plt.show()
 
很简单，代码分为两部分，上面为主要的计算系数w的值，下面为直观查看每个次项拟合后的效果图，如下：
 
 
可以清楚的看到当3次项（也就是m=3）时，效果还是一般，但从4次项（m=4）后，对于数据的拟合性就明显优于3次项了，当 m=8 时，曲线呈现出明显的震荡，这也就是线性回归实验中所讲到的过拟和（Overfitting）现象，后面再介绍如何解决这个问题。
 
使用scikit-learn进行多项式拟合：
 
对于一个二次多项式而言，我们知道它的标准形式为：，但是，多项式回归其实相当于是线性回归的特殊形式（开头也提到）。例如，我们这里令，，那么原方程就转换为：，这也就变成了多元线性回归。完成了一元高次多项式到多元一次项式之间的转换。（如下，看作将多元一次项合并为一个矩阵中线性求解）
 
举例说明，对于自变量向量  和因变量 ，如果 ：
 
                      [[ 3]
                  [-2]
                        [ 4]]
 
则可以通过线性回归模型进行拟合。同样，如果对于一元二次多项式，能得到，构成的特征矩阵，即：
 
                                        [[ 3.  9.]
                 [-2.  4.]
                                          [ 4. 16.]]
 
那么也就可以使用线性回归进行拟合了。
 
所以这里有了一个方便的工具，scikit-learn 中，我们可以通过 PolynomialFeatures() 类自动产生多项式特征矩阵
 
 
 
sklearn.preprocessing.PolynomialFeatures(degree=2, interaction_only=False, include_bias=True)
 
 
 
degree: 多项式次数，默认为 2 次多项式
 
 
interaction_only: 默认为 False，如果为 True 则产生相互影响的特征集。
 
 
include_bias: 默认为 True，包含多项式中的截距项。
 
通过下面代码解决上面一样的问题：
 
"""
	使用 PolynomialFeatures 自动生成特征矩阵
"""
from sklearn.preprocessing import PolynomialFeatures

x = np.array(x).reshape(len(x), 1) # 转换为列向量
y = np.array(y).reshape(len(y), 1)

poly_features = PolynomialFeatures(degree=2, include_bias=False)    # 特征矩阵模型构建用来多项式的特征矩阵
poly_x = poly_features.fit_transform(x)    # 通过模型转换x数据

"""
    用线性回归拟合上面转换后的数据
"""
from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(poly_x,y)		# 这里使用转换后的特征矩阵相当于将一元2次转换为二元一次
# print(model.intercept_,model.coef_)

"""绘制拟合图像
"""
x_temp = np.linspace(0, 80, 10000)
x_temp = np.array(x_temp).reshape(len(x_temp),1)
poly_x_temp = poly_features.fit_transform(x_temp)

plt.plot(x_temp, model.predict(poly_x_temp), 'r')
plt.scatter(x, y)
plt.show()         # 得到的结果与上面用leastsq的二次项结果一致
 
 
这里二次项的图（degree=2），若设置更高的次项3，4，5（degree=3，4，5），则可以得到上面leastsq方法相同的效果，只需要通过PolynomialFeatures方法生成对应次项的特征矩阵就行。
 
评价指标
 
可以使用前一篇文章提到的 平均绝对误差（MAE）、均方误差（MSE）等方法来衡量，具体可以使用sklearn中的以下两个方法：
 
 
 
from sklearn.metrics import mean_absolute_error
 from sklearn.metrics import mean_squared_error
 
 
2次多项式不一定比线性回归的效果好，需要根据具体的测试结果做决定，可以试试更高次的多项式回归预测 ，下面我们便来验证以一下。
 
make_pipeline通道：
 
通过实例化 make_pipeline 管道类，实现调用一次 fit 和 predict 方法即可应用于所有预测器。make_pipeline 是使用 sklearn 过程中的技巧创新，其可以将一个处理流程封装起来使用。 更详细的使用可以参考这里
 
上面的多项式回归中，我们需要先使用 PolynomialFeatures 完成特征矩阵转换，再放入 LinearRegression 中。那么，PolynomialFeatures + LinearRegression 这一个处理流程，就可以通过 make_pipeline 封装起来使用。
 
以下数据集并非上面的x，y，如下数据（将Year作为特征列train_x，Values为目标值train_y，且数据是按0.7比例分割为训练、测试集）：
 
 
"""更高次多项式回归预测
"""
from sklearn.pipeline import make_pipeline

train_x = train_x.reshape(len(train_x),1)
test_x = test_x.reshape(len(test_x),1)
train_y = train_y.reshape(len(train_y),1)

for m in [3, 4, 5]:
    model = make_pipeline(PolynomialFeatures(m, include_bias=False), LinearRegression())
    model.fit(train_x, train_y)
    pre_y = model.predict(test_x)
    print("{} 次多项式回归平均绝对误差: ".format(m), mean_absolute_error(test_y, pre_y.flatten()))
    print("{} 次多项式均方根误差: ".format(m), mean_squared_error(test_y, pre_y.flatten()))
    print("---")
 
 
可以看到高次项的误差都比线性回归要低，而2次项则要高，所以对于模型均需要进行验证后使用。
 
好了，上面介绍了 多项式回归的原理与Scipy实现、sklearn构建特征矩阵转为线性运算、Pipeline方便的通道运算、以及误差评判，但还有一个很重要的问题，那就是如何选择最佳的n次项？
 
几次项才是最佳的选择？
 
其实这个问题很简单，我们可以设置一个误差指标（MSE、MAE等），然后绘制增加次项后，预测结果的误差值图。观察图形来选择一个合适的次项点（类似肘部法则也可以作为参考）：
 
"""
	计算 m 次多项式回归预测结果的 MSE 评价指标并绘图
"""
mse = [] # 用于存储各最高次多项式 MSE 值
m = 1 # 初始 m 值
m_max = 10 # 设定最高次数
while m <= m_max:
    model = make_pipeline(PolynomialFeatures(m, include_bias=False), LinearRegression())
    model.fit(train_x, train_y) # 训练模型
    pre_y = model.predict(test_x) # 测试模型
    mse.append(mean_squared_error(test_y, pre_y.flatten())) # 计算 MSE
    m = m + 1

# print("MSE 计算结果: ", mse)
# 绘图
plt.plot([i for i in range(1, m_max + 1)], mse, 'b')
plt.scatter([i for i in range(1, m_max + 1)], mse)

# 绘制图名称等
plt.title("MSE of m degree of polynomial regression")
plt.xlabel("m")
plt.ylabel("MSE")
plt.show()
 
 
如上图，可以明显看到当次数为3时，误差指标已经趋近于平稳了，所以这里选择3最好，当次数项越多可能会出现过拟合的问题，模型的泛化能力会降低。
 
 
 
参考文章：
 
https://baike.baidu.com/item/%E5%A4%9A%E9%A1%B9%E5%BC%8F%E5%9B%9E%E5%BD%92/21505384?fr=aladdin
 
https://blog.csdn.net/lanchunhui/article/details/50521648

多项式的各种操作
 
By SemiWaker
 
目录
 
多项式乘法点这里
多项式逆元
多项式除法
多项式牛顿迭代法
多项式对数
多项式指数函数
多项式幂函数
多项式三角函数
多项式多点求值
多项式快速插值 
 
  
 
多项式逆元
 
定义
 

 $$A(x)A^{-1}(x)\equiv 1\pmod {x^n}$$ 
 
A(x)：要求逆元的多项式
$A^{-1}(x)$：多项式 A(x) 的逆元
$\pmod{x^n}$：在模$x^n$的意义下。$x^n$是一个多项式，不是一个数。 
 模$x^n$的意义类似数的除法：设$a=bd+r$，且$r<b$，那么有$a\equiv r\pmod b$。在多项式中，只要r的界比b的界小就好了。 
 多项式除法要用多项式逆元，但是多项式逆元的定义又是多项式除法，这里先介绍多项式逆元，所以除法部分简单略过 
 $x^n$比较特殊，因为模$x^n$实际上就是舍弃掉$\ge n$次的所有项，只剩下 0 到 n-1次的项。 
 （原因显然$\sum_{i=0}^{N-1} a_ix^i=x^n\sum_{i=0}^{N-n-1} a_ix^i+\sum_{i=0}^{n-1} a_ix^i$） 
 为什么要模$x^n$呢？ 
 因为如果不模的话，除了只有常数项的多项式以外，其他多项式的逆元都是无穷级数。 
 例：$1+2x+3x^2$。为了方便我们只写系数数列，从次数小的开始。 
 可以写成 1 1 1 。 
 设这个数列是$a_i$，然后设逆元的系数数列是$b_i$，最后设它们乘起来是$c_i$。 
 根据定义，$c_i$满足$c_0=1$，$\forall i>0,c_i=0$。 
 然后我们考虑这个的逆元会是什么。 
 显然$b_0=1$。 
 然后$c_1=a_0b_1+a_1b_0$，可以解出$b_1=-1$ 
 然后$c_2=a_0b_2+a_1b_1+a_2b_0$，可以解出$b_2=0$ 
 依次类推$b_3=1$，$b_4=-1$…… 
 发现无论如何都有非0项存在，所以是无穷级数。 
 模$x^n$之后就可以解决这个问题，我们就可以只求我们需要的前n项。
 
做法
 
暴力做法
 
上面举例为什么一定要模$x^n$的时候，已经给出了暴力做法了。 
 直接递推，时间复杂度$O(n^2)$
 
非暴力：
 
基本思路：倍增 
 只有常数项的多项式求逆元是比较简单的事，即普通的求逆元。 
 然后考虑怎么让已经求出的多项式翻倍，也就是倍增的思路。
 
首先，我们设长度为$\lceil\frac n2\rceil$的已经求出的逆元为B(x)。 
 注意，是向上取整。 
 向上取整的好处是：翻倍后保证长度$\ge n$。
 

 $$A(x)B(x)\equiv 1\pmod {x^{\lceil\frac n2\rceil}}$$  

 然后直接设当前求的逆元为

$A^{-1}(x)$ 

 

 $$A(x)A^{-1}(x)\equiv 1\pmod{x^n}$$ 
 
注意到，前n项为(1 0 0 0 0 ……），那么前$\lceil\frac n2\rceil$项也是(1 0 0 0 ……）。 
 所以当前的逆元折半后还是逆元。 
 
 $$A(x)A^{-1}(x)\equiv 1\pmod {x^{\lceil\frac n2\rceil}}$$  

 把两个式子合在一起得到

 $$A(x)(B(x)-A^{-1}(x))\equiv 0\pmod{x^{\lceil\frac n2\rceil}}$$  

 显然A(x)不为0，所以后面的为0。 

 

 $$B(x)-A^{-1}(x)\equiv 0\pmod{x^{\lceil\frac n2\rceil}}$$ 
 
我们现在有了之前的逆元和当前的逆元的关系，但是这时在模$\lceil\frac n2\rceil$的情况下的，考虑怎样才能变成在模n的情况下表示之前逆元和当前逆元的关系。
 
有一个很简单的方法：把它平方了。 
 
 $$(B(x)-A^{-1}(x))^2\equiv 0\pmod{x^n}$$  

 为什么一个在模

$\lceil\frac n2\rceil$的情况下为0的多项式，平方之后模

$x^n$也是0呢？ 

 因为模

$\lceil\frac n2\rceil$的情况下为0的意思是前

$\lceil\frac n2\rceil$项都是0。 

 平方之后，显然前

$\lceil\frac n2\rceil$项也还是0。 

 然后我们考虑后面的那些项。 

 因为是卷积，所以

$a_n=\sum a_ia_{n-i}$。其中 i 和 n-i 必然有一个小于等于

$\lceil\frac n2\rceil$，那么每一项至少有一个数是0，所以全部都是0。 

 那么就说明，平方之后的多项式的前n项必然都是0，所以我们可以用平方来倍增。
 
接着把平方拆开。 
 
 $$B^2(x)+A^{-2}(x)-2B(x)A^{-1}(x)\equiv 0\pmod {x^n}$$  

 

$A^{-2}(x)$会很不好处理，我们整体乘一个A(x)。（显然A(x)非0，所以不会有问题。） 

 得到： 

 

 $$B^2(x)A(x)+A^{-1}(x)-2B(x)\equiv 0\pmod {x^n}$$  

 这样我们就凑出了现在的逆元和之前的逆元的关系了。 

 

 $$A^{-1}(x)\equiv B(x)(2-B(x)*A(x))\pmod {x^n}$$  

 注意到B(x)只有

$\lceil\frac n2\rceil$这么多项，所以后面的式子应该最多2n项。 

 所以把B(x)和A(x)用FFT得到2n项的点值然后计算即可。
 
时间复杂度 
 
 $$T(n)=T(\frac n2)+O(nlogn)=O(nlogn)$$  

 常数大概是FFT的6倍。（每一层3个FFT，加起来要遍历两次。）
 
代码
 
void PolyRev(int n,Complex *A,Complex *B)
{
    static Complex tmp[MAXN*4];
    if (n==1) {B[0]=1.0/A[0];return;}
    PolyRev((n+1)>>1,A,B);
    int N=1;while (N<n<<1) N<<=1;
    for (int i=0;i<n;++i) tmp[i]=A[i];
    for (int i=n;i<N;++i) tmp[i]=0;
    FFT(N,tmp,1);FFT(N,B,1);
    for (int i=0;i<N;++i) B[i]=B[i]*(2.0-tmp[i]*B[i]);
    FFT(N,B,-1);
    for (int i=n;i<N;++i) B[i]=0;
}
 
NTT版类似，这里用复数只是方便理解
 
注意点
 
是的，n不一定是2的幂，因为每次向上取整，所以n是多少都可以。
多项式是否有逆元，根据上面的过程，只取决于常数项是否有逆元。 
 
  
 
多项式除法
 
定义
 

 $$A(x)=B(x)D(x)+R(x)$$ 
 
A(x)：被除数多项式，界为n。
B(x)：除数多项式，界为m，显然有$m<n$。
D(x)：商多项式，界为n-m+1。
R(x)：余数多项式，界小于m。 
 类似普通的除法中余数要小于除数才能保证结果唯一，虽然多项式无法比较大小，但是R(x)的界比B(x)的界小，就可以保证结果唯一了。
 
多项式取模
 

 $$A(x)\equiv R(x)\pmod{B(x)}$$  

 即上面除法式子的模的形式，是等价的。
 
做法
 
思路：转换 
 我们观察除法和逆元的差别： 
 1、有余数，所以$B(x)D^{-1}(x)\neq 1\pmod {x^n}$ 
 2、没有模，不需要模除法的结果长度也是确定的。 
 我们知道逆元怎么做，那么我们试着用逆元来做除法。
 
先考虑求出商，求出商之后余数是可以带进去求出来的。 
 第二个差别很容易解决，我们只需要强行模$x^{n-m+1}$就可以保证除法的结果不变。 
 那么有余数的问题怎么办呢？ 
 能不能在模的时候把余数去掉？ 
 观察每一个多项式的项数： 
 A(x) n 
 B(x) m 
 D(x) n-m+1 
 R(x) <=m-1 
 也就是说，在模$x^{n-m+1}$的情况下，会抛弃掉前m-1项。 
 但是，R(x) 只有后m-1项，所以直接模无法去掉余数。
 
考虑怎么让后m-1项移动到前m-1项。 
 定义反转操作： 
 
 $$A^R(x)=x^{n-1}A(\frac 1x)=\sum_{i=0}^na_{n-i-1}x^i$$  

 也就是，我们把 

$\frac1x$ 带入 A(x)，然后乘上

$x^{n-1}$。 

 这个操作的意义是什么呢？ 

 
反转系数。 

 原来的某一项

$a_ix^i$，变为

$a_ix^{-i}\times x^{n-1}=a_ix^{n-i-1}$。 

 也就是系数反转了。
 
我们试着用系数反转来去掉余数。 
 
 $$A^R(x)=x^{n-1}A(\frac1x)=x^{m-1}B(\frac1x)x^{n-m}D(\frac1x)+x^{n-m+1}x^{m-2}R(\frac1x)$$  

 也就是： 

 

 $$A^R(x)=B^R(x)D^R(x)+x^{n-m+1}R^R(x)$$  

 这时再模

$x^{n-m+1}$ 

 

 $$A^R(x)\equiv B^R(x)D^R(x)\pmod{x^{n-m+1}}$$ 
 
然后问题解决。 
 具体步骤如下： 
 反转A(x)、B(x)，求$B^R(x)$在模$x^{n-m+1}$意义下的逆元，然后乘上$A^R(x)$得到$D^R(x)$，再反转过来就得到了商D(x)。 
 如果要求R(x)，代入$R(x)=A(x)-B(x)D(x)$即可。
 
需要的只有逆元和乘法而已。 
 时间复杂度O(nlogn)
 
注意点
 
1、除法没有模，不要写多了。 
 2、要保证$m\le n$，否则没有意义。 
 3、所有计算长度不会超过2n，所以FFT的长度是2n。
 
程序
 
void PolyDiv(int n,Complex *A,int m,Complex *B,Complex *D,Complex *R)
{
    static Complex A0[MAXN*4],B0[MAXN*4],B1[MAXN*4];
//D
    int N=1;while (N<n<<1) N<<=1;
    for (int i=0;i<N;++i) A0[i]=B0[i]=B1[i]=D[i]=R[i]=0.0;
    for (int i=0;i<n;++i) A0[n-i-1]=A[i];
    for (int i=0;i<m;++i) B0[m-i-1]=B[i];
    PolyRev(n-m+1,B0,B1);
    FFT(N,A0,1);FFT(N,B1,1);
    for (int i=0;i<N;++i) D[i]=A0[i]*B1[i];
    FFT(N,D,-1);
    for (int i=n-m+1;i<N;++i) D[i]=0.0;
    for (int i=0,j=n-m;i<j;++i,--j) swap(D[i],D[j]);
//R
    for (int i=0;i<m;++i) B0[i]=B[i];
    for (int i=0;i<n-m+1;++i) B1[i]=D[i];
    for (int i=n-m+1;i<N;++i) B1[i]=0.0;
    FFT(N,B0,1);FFT(N,B1,1);
    for (int i=0;i<N;++i) R[i]=B0[i]*B1[i];
    FFT(N,R,-1);
    for (int i=0;i<n;++i) R[i]=A[i]-R[i];
}
 
 
多项式牛顿迭代法
 
定义
 
普通的牛顿迭代法
 
问题：求函数f(x)的零点。 
 首先随便选一个点x0。 
 然后把f(x)泰勒展开。 
 
 $$f(x)=f(x_0)+f'(x_0)(x-x_0)+\frac {f''(x_0)(x-x_0)^2}2+\cdots$$  

 为了简单，只保留线性部分。 

 

 $$f(x)=f(x_0)+f'(x_0)(x-x0)=0$$  

 然后得到 

 

 $$x=x0-\frac{f(x_0)}{f'(x_0)}$$  

 当然，这个不是准确解。 

 只需要用x代替x0，迭代多次后就可以得到近似解。
 
多项式牛顿迭代法
 
问题：求函数g(f(x))的零点多项式法f(x) 
 解释：g(f(x))是一个关于多项式f(x)的函数，其中f(x)是变量。 
 即求一个多项式f(x)使得
 $$g(f(x))\equiv0\pmod{x^n}$$ 
 
举例
 

 $$g(f(x))=3f^2(x)-2f(x)+1$$  

 这个函数的零点即方程

$3f^2(x)-2f(x)+1=0$的解， 

 

 $$g(f(x))=f(x)-A^k(x)$$  

 这个函数的零点即

$A^k(x)$，也就是把A(x)的k次方求出来。 

 

 $$g(f(x))=lnf(x)-A(x)$$  

 这个函数的零点为

$e^{A(x)}$，即A(x)的exp。 

 
简单来说，多项式牛顿迭代法就是用来解关于多项式的方程的。
 
做法
 
思路：倍增 
 当n=1的时候，只有常数项，可以直接算出。 
 常数项的解可能有多个，代入哪一个决定了最后算出哪个多项式 
 然后考虑长度$\lceil\frac n2\rceil$的多项式$f_0(x)$和长度为n的$f(x)$的关系。 
 我们可以把$f_0(x)$当成普通的牛顿迭代法中的$x_0$去泰勒展开，然后观察会怎样： 
 
 $$g(f(x))=g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))+\frac{g''(f_0(x))}2(f(x)-f_0(x))^2+\cdots$$  

 注意到：

$f_0(x)$和

$f(x)$的前

$\lceil\frac n2\rceil$项必然是相同的。 

 那么，

$f(x)-f_0(x)$的前

$\lceil\frac n2\rceil$项也会是0。 

 所以

$(f(x)-f_0(x))^2$以及更高次方的前n项都是0。 

 所以，
对于多项式来说，只保留线性部分得到的是准确解。 

 然后我们列出对于多项式的牛顿迭代公式： 

 

 $$f(x)=f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}$$  

 代入求值即可。
 
每一层需要逆元和乘法，所以时间复杂度: 
 
 $$T(n)=T(\frac n2)+O(nlogn)=O(nlogn)$$  

 常数的话，每一层需要至少3次FFT和一次逆元，相当于9次FFT，所以总共大约18次FFT，即常数至少是FFT的18倍。 

 
几乎等于多了一个logn
 
注意点
 
1、多项式函数求导的时候，常数多项式应该视为常数去求导，不应该把常数多项式直接求导。 
 例如$g(f(x))=f^2(x)-A(x)$ 
 求导得：$g'(f(x))=2f(x)$ 
 而不是$2f(x)-A'(x)$ 
 2、FFT的长度是2n。 
 3、精准度取决于常数项，之后的操作都是完全没有误差的。（如果无视FFT的话。） 
 4、多个零点的函数，求出哪个零点取决于常数项代入哪个零点。零点的数量也是取决于常数项。
 
 
多项式对数
 
定义
 
对于一个多项式A(x)，求： 
 
 $$lnA(x)\pmod{x^n}$$  

 即一个n项的多项式。 

 
说实话，多项式对数到底是什么我并不清楚，但是知道功能和对数一样就可以了。
 
做法
 
直接通过数学的方法算出。 
 即 
 
 $$lnA(x)=\int (lnA(x))'=\int\frac{A'(x)}{A(x)}$$  

 
为啥这里又直接把A(x)求导了呢？因为这里的A(x)不是常数。 

 然后步骤就很简单了： 

 求出A(x)在模

$x^n$意义下的逆元，再求出

$A'(x)$，乘起来之后再积分。
 
多项式的求导和积分都是可以O(n)完成的。 
 简单来说： 
 假设A(x)的系数数列为ai，求导之后的系数数列为bi，积分之后的为ci。 
 则有： 
 
 $$b_i=(i+1)a_{i+1}$$  

 

 $$c_i=\frac {a_{i-1}}i$$  

 记得将空出来的一项设为0就好了。
 
时间复杂度O(nlogn)
 
程序
 
void PolyLn(int n,Complex *A,Complex *B)
{
    static Complex A0[MAXN*4],A1[MAXN*4];
    int N=1;while (N<n<<1) N<<=1;
    for (int i=0;i<N;++i) A0[i]=0.0,A1[i]=0.0;
    PolyRev(n,A,A0);
    for (int i=0;i<n-2;++i) A1[i]=A[i+1]*double(i+1);A1[n-1]=0.0;
    FFT(N,A0,1);FFT(N,A1,1);
    for (int i=0;i<N;++i) A0[i]*=A1[i];
    FFT(N,A0,-1);
    for (int i=0;i<n-2;++i) B[i+1]=A0[i]/double(i+1);B[0]=0.0;
    for (int i=n;i<N;++i) B[i]=0.0;
}
 
 
多项式指数函数
 
定义
 
对于一个多项式A(x)，求： 
 
 $$e^{A(x)}\pmod{x^n}$$  

 即exp(A(x))。
 
做法
 
牛顿迭代 
 首先数学变换 
 
 $$f(x)\equiv e^{A(x)}\pmod{x^n}$$  

 两边取ln。 

 

 $$lnf(x)-A(x)\equiv0\pmod{x^n}$$  

 定义多项式函数g(f(x))为： 

 

 $$g(f(x))\equiv lnf(x)-A(x)\pmod{x^n}$$  

 然后带入牛顿迭代公式 

 

 $$f(x)=f_0(x)-\frac{lnf_0(x)-A(x)}{\frac 1{f_0(x)}}$$  

 化简得： 

 

 $$f(x)=f_0(x)(1-lnf_0(x)+A(x))$$ 
 
然后按照牛顿迭代的方法去做就好了。 
 时间复杂度O(nlogn)
 
程序
 
void PolyExp(int n,Complex *A,Complex *B)
{
    if (n==1) {B[0]=exp(A[0]);return;}
    static Complex A0[MAXN*4],B0[MAXN*4];
    PolyExp((n+1)>>1,A,B);
    int N=1;while (N<n<<1) N<<=1;
    for (int i=0;i<N;++i) B0[i]=0.0;
    PolyLn(n,B,B0);
    for (int i=0;i<n;++i) A0[i]=A[i];for (int i=n;i<N;++i) A0[i]=0.0;
    FFT(N,A0,1);FFT(N,B,1);FFT(N,B0,1);
    for (int i=0;i<N;++i) B[i]=B[i]*(1.0-B0[i]+A0[i]);
    FFT(N,B,-1);
    for (int i=n;i<N;++i) B[i]=0.0;
}
 
 
多项式幂函数
 
定义
 
对于一个多项式A(x)和一个数k，求： 
 
 $$A^k(x)\pmod{x^n}$$ 
 
做法
 
数学变换 
 
 $$A^k(x)\equiv e^{klnA(x)}\pmod{x^n}$$  

 只需要先求

$lnA(x)$，再乘上k，再求exp即可。 

 实际，乘上一个多项式也是可行的。所以用这种方法也可以求出多项式的多项式次方。 

 时间复杂度O(nlogn)。 

 
当然，常数太大，如果是正整数的话直接快速幂也差不多了。
 
程序
 
void PolyPow(int n,Complex *A,Complex *B,double k)
{
    static Complex A0[MAXN*4];
    PolyLn(n,A,A0);
    for (int i=0;i<n;++i) A0[i]*=k;
    PolyExp(n,A0,B);
}
 
 
多项式三角函数
 
定义
 
对于一个多项式A(x)，求： 
 
 $$sin(A(x))\pmod{x^n}$$  

 

 $$cos(A(x))\pmod{x^n}$$  

 求出这两个三角函数之后就可以求出所有三角函数了。 

 
虽然意义不明
 
做法
 
求导是没有什么用的，求导之后还是三角函数。 
 所以需要找一个三角函数和其他初等函数的关系式。 
 欧拉公式 
 
 $$e^{ix}=cos(x)+isin(x)$$  

 把A(x)带入x得到 

 

 $$e^{iA(x)}=cos(A(x))+isin(A(x))$$  

 这样一次exp就可以求出cos(A(x))和sin(A(x))。 

 唯一需要注意的问题就是乘了i之后怎么做exp。 

 首先，多项式的系数要用复数表示。 

 然后考虑FFT，实际上FFT的过程完全没有变化，直接代入即可。 

 接着考虑逆元，没有什么区别，常数项的逆元也是

$\frac 1{a_0}$，只不过是复数的除法而已。 

 然后考虑ln，完全没区别。 

 最后考虑exp，有区别的地方在于常数项，要代入普通的欧拉公式来求。 

 这时求出来的显然是无理数，所以NTT就不要想了。
 
 
多项式多点求值
 
多项式快速插值
 
留坑待填
 
 
BySemiWaker

 一般情况下，一元n次多项式可写成：
 
 
 
 
 其中，pi是指数为ei的项的非零系数，且满足
 
 
 
 
 
 因此，我们可以采用线性表（定义：线性表是由n个数据元素构成的有限序列，比如数组、向量、链表等等）来表示：
 
 
 
 
 
 
 
 其中，每一项的指数i可以用其系数pi的序号表示。
 
        在通常的应用中，多项式的次数比较大，使得线性表的长度很难确定，因此我们可以考虑链表，向量也可以（c++中）。举例说明：假如我们用数组来表示下面的多项式：
 


 
        可见，我们需要一个大小为1549的数组来表示，而实际有用的信息只有数组中的4个元素，其他地方都是0，所以造成了空间浪费。并且如果我们事先不知道多项式的最高次项的指数，则我们需要定义一个足够大的数组来存储，这样做显然浪费了很多内存空间。我们可以使用链表来解决上述问题：
 
 
 
        在计算机内，我们用一个结点来存放多项式的一项，为了节约空间，并和书写习惯一致，只需保留非零系数的项。每个结点分系数、指数和指针三个域，如下图所示，其中的指针next指明下一项的位置。
 
 
 
 
 
 
 例如，下面多项式分别为A，B:
 
         
 
 用循环链表可以表示如下：
 
 
 
 
        两个多项式相加的运算规则很简单，对所有指数相同的项，将其对应系数相加，若和不为零，则构成和多项式中的一项；将所有指数不相同的项复制到和多项式中。具体实现时，我们以上面的多项式A，B为测试样例。可采用另建链表来存储和的多项式的方法，或采用把一个多项式归并入另一个多项式的方法。我们以后种方法为例,即将A+B的和多项式存储到A中。具体程序实现如下（我采用了循环链表）：
 
 
 
 
 
 
 
  
 
   [cpp] 
   view plain
    copy 
   
 
   
 
    
   
   
  
 
 
 
 
 #include<stdio.h>  
 #include<stdlib.h>  
 #include<conio.h>  
   
 typedef struct pnode//用链表来存储多项式信息  
 {  
     float coef;//多项式系数  
     int   exp;//多项式指数  
     struct pnode *next;  
 }polynode;  
   
 polynode *Create()  
 {  
     float coef;  
     int exp;  
     polynode *head,*s,*r;  
     head=(polynode*)malloc(sizeof(polynode));  
     head->coef=0;  
     head->exp=-1;  
     r=head;  
     printf("请输入各项的系数和指数：\n");  
     while(1)  
     {  
         scanf("%f %d",&coef,&exp);  
         if(coef!=0)//输入0 0来结束输入  
         {  
             s=(polynode*)malloc(sizeof(polynode));  
             s->coef=coef;//s用来保存当前节点  
             s->exp=exp;  
             r->next=s;  
             r=s;  
         }  
         else  
             break;  
     }  
   
     r->next=head;//构造循环链表  
     return head;  
   
 }  
   
 polynode*PolyAdd(polynode* pa,polynode* pb)//进行多项式相加  
 {  
     polynode *p,*q,*r,*s;  
     float x;  
     p=pa->next;//分别指向多项式的第一项  
     q=pb->next;  
     s=pa;//s用于保存当前节点  
   
     while((p!=pa)&&(q!=pb))//没有结束，回到链表头  
     {  
         if(p->exp<q->exp)//p的指数小于q的指数，将p放入链表中  
         {  
             s=p;  
             p=p->next;  
         }  
         else if(p->exp>q->exp)//p的指数大于q的指数，将q放入链表中  
         {  
             r=q->next;  
             q->next=p;  
             s->next=q;  
             s=q;  
             q=r;  
         }  
         else//当两者指数相同时，进行合并  
         {  
             x=p->coef+q->coef;  
             if(x!=0)  
             {  
                 p->coef=x;  
                 s=p;  
             }  
             else//若合并结果为0，将该节点移除  
             {  
                 s->next=p->next;  
                 free(p);  
             }  
                     p=s->next;  
         r=q;  
         q=q->next;  
         free(r);  
         }  
   
     }  
   
     if(q!=pb)//如果多项式b的项数少于多项式a的情况  
     {  
         r=q;  
         while(r->next!=pb)  
             r=r->next;  
         s->next=q;  
         r->next=pa;  
     }  
     return pa;  
 }  
   
 void Output(polynode *head)// 输出多项式信息  
 {  
     polynode *p;  
     printf("系数和指数分别为：");  
     p=head->next;  
     while(p!=head)  
     {  
         printf("%.1f , %d    ",p->coef,p->exp);  
         p=p->next;  
     }  
     printf("\n");  
 }  
   
 void main()  
 {  
     polynode* ha,*hb;  
     printf("\n建立多项式A:");  
     ha=Create();  
     Output(ha);  
     printf("\n建立多项式B:");  
     hb=Create();  
     Output(hb);  
   
     ha=PolyAdd(ha,hb);  
     printf("\n多项式A+B：");  
     Output(ha);  
 }  
 
 
 运行结果如下：
 
 
 
 
 
 
 
 原文：http://blog.csdn.net/tengweitw/article/details/40476935
 
 作者：nineheadedbird

多项式乘法入门
 
By SemiWaker
 
这是一篇蒟蒻对FFT、DFT、CZT、NTT的弱鸡理解
 
多项式
 

 $$a_0x+a_1x^1+a_2x^2\cdots+a_{n-1}x^{n-1}$$  

 上面的这个形式叫做多项式。 

 系数：

$a_{0..n-1}$ 

 项：

$a_ix^i$ 

 界：n 

 为了方便我们系数序列就可以表示多项式。
 
线性卷积
 

 $$A\times B=\sum_{i=0}^{2n-2}(\sum_{j=0}^{i}A_jB_{i-j})x^i$$  

 简单来说，就是把两个多项式直接乘起来。 

 第i项的构成如下： 

 任取两个项

$A_jx^j$和

$B_kx^k$，如果

$j+k=i$，那么得到

$A_jB_kx^i$。 

 将所有

$j+k=i$的项系数两两相乘加起来即可。
 
循环卷积
 

 $$A\times B=\sum_{i=0}^{n-1}(\sum_{j,k}A_jB_k[j+k\equiv i\pmod n])x^i$$  

 比较难理解。 

 其实就是把多项式的-1项设为n-1项，-2项设为n-2项…… 

 然后，线性卷积为了防止出现-1项，设定了j<=i，现在我们把它去掉。 

 

 $$A\times B=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}A_jB_{i-j})x^i$$ 
 
线性卷积和循环卷积关系
 
我们用比较形象地方式说明。 
 举个例子 
 
 $$(1+2x+3x^2+4x^3)\times(5+6x+7x^2+8x^4)$$  

 即

$(1,2,3,4)\times(5,6,7,8)$ 

 线性卷积： 

 
 

 循环卷积 

 
 

 如果把 (5,6,7,8) 复制展开 

 
 

 然后我们观察，设线性卷积的结果为

$C_{0...6}$，循环卷积的结果为

$D_{0..3}$ 

 则有： 

 

$D_0=C_0+C_4$ 

 

$D_1=C_1+C_5$ 

 

$D_2=C_2+C_6$ 

 

$D_3=C_3$ 

 也就是说，把线性卷积重叠在一起，就得到了循环卷积。
 
多项式点值表示
 
我们把N个数$x_{0..n-1}$带入x，可以得到多项式的N个值$A(x_{0..n-1})$。 
 将每一个数和对应的值当成一个点$(x_i,A(x_i))$，这N个点叫做多项式的点值表示。 
 为什么点值表示可以代表多项式呢? 
 因为我们可以反过来用点值表示求出多项式。 
 待定系数+解方程就可以了。
 
多项式插值
 
就是把点值表示转换成系数表示的一个过程。
 
DFT
 
离散傅里叶变换 
 考虑怎样快速求多项式卷积。 
 按照定义去写$O(n^2)$，显然不优。
 
我们可以考虑点值表示。 
 两个多项式如果取值用的是同样的数$x_{0..n-1}$，那么得到的值直接乘起来就可以得到卷积之后的多项式的点值表示。
 
但是如果按照定义一个一个数带入求多项式的值，还是$O(n^2)$的，没有优化。
 
我们考虑带入特殊的数去求点值表示。
 
DFT就是这样一个过程：将一个多项式转化为用n次单位根表示的点值表示。
 
FFT是实现DFT的算法。
 
单位根
 
简单来说$x^n=1$的复数解。 
 n次单位根有n个复数解，设为$\omega_n^k$，其中k=0…n-1。 
 $\omega_n^k=e^{\frac {2k\pi I}n}$ 
 用欧拉公式展开： 
 $e^{xI}=cos(x)+sin(x)I$ 
 得到 
 $(\omega_n^k)^n=cos(2k\pi)+sin(2k\pi)I=1$
 
画在复数平面上，刚好是把单位元n等分的n个点。
 
有一些有趣的性质 
 $\omega_n^k=\omega_n^{k\pmod n}$ 
 由三角函数周期性可证。 
 $\omega_n^{\frac n2}=-1$ 
 显然。 
 $\omega_n^n=1$ 
 显然。 
 $\omega_n^k=\omega_{\frac n2}^{\frac k2}$ 
 显然。 
 $(\omega_n^k)^2=(\omega_n^{k+\frac n2})^2$ 
 把定义带入可证。
 
FFT的实现
 
快速傅里叶变换 Cooley-Tukey算法 
 考虑把$\omega_n^{0...n-1}$一起带入求值。
 
则
 $$A_k=\sum_{j=0}^{n-1}a_j(\omega_n^k)^j$$  

 可以分治为两部分 

 

 $$A_k=\sum_{j=0}^{\frac {n-1}2}a_{2j}(\omega_n^k)^{2j}+(\omega_n^k)\sum_{j=0}^{\frac {n-1}2}a_{2j+1}(\omega_n^k)^{2j}$$  

 两半的长度都为

$\frac n2$
 
而$(\omega_n^k)^{2j}=(\omega_n^{k+{n/2}})^{2j}=\omega_{\frac n2}^{kj}$ 
 所以而$\omega_{\frac n2}^{kj}$只有$\frac n2$个取值，而且每个出现两次。
 
这样就通过分治减小了规模。
 
现在变成了对两个多项式$(a0,a2,a4...)$和$(a1,a3,a5...)$求带入$\omega_{\frac n2}^{0..\frac n2-1}时的值$。
 
分治完了之后，我们要求回原来的多项式。 
 设分治的结果为$B_{0..\frac n2-1}$和$C_{0..\frac n2-1}$。 
 则$A_k=B_k+\omega_n^kC_k$ 
 但是注意此时k只有$\frac n2$ 
 所以$A_{k+\frac n2}=B_k+\omega_n^{k+\frac n2}C_k$ 
 又$\omega_n^{k+\frac n2}=-\omega_n^k$ 
 所以$A_{k+\frac n2}=B_k-\omega_n^kC_k$
 
用一个简单的图来记： 
  
 设当前分治长度为L，则 
 $A_k=A_k+\omega_n^kA_{k+\frac L2}$ 
 $A_{k+\frac L2}=A_k-\omega_n^kA_{k+\frac L2}$ 
 注意变量自我迭代时，要开临时变量 
 这两个位置刚好交错计算，形状类似蝴蝶，所以叫做蝴蝶变换。
 
边界： 
 n=1时，$\omega_1^1=1$，所以直接把$a_k$放进去就好了。
 
程序怎么实现呢？ 
 直接分治是可以的，但是有更好的方法。 
 我们考虑分治时的分类方法： 
 每一层按照奇偶数。 
 如果我们保持编号不变，那么就变成：第i行按照从低到高第i位分。 
 由于是从低到高，所以最后的排列为：每一个数的位置，为编号二进制倒序后位置。 
 举例：有8个数，编号0~7。 
 分治过程： 
 000 001 010 011 100 101 110 111 
 000 010 100 110|001 011 101 111 
 000 100|010 110|001 101|011 111 
 000|100|010|110|001|101|011|111
 
把最后一行的二进制倒过来： 
 000 001 010 011 100 101 110 111 
 刚好是从0~7。
 
设Reverse(x)为x二进制倒序后的数 
 一开始我们可以将ai放到Reverse(i)的位置。
 
然后设当前层分治长度为L。 
 每L个一起处理，进行蝴蝶变换即可。
 
进一步优化，注意蝴蝶变换中$\omega_n^k$的k取值为0~L/2 
 为了尽量减少对$\omega_n^k$的计算次数，我们可以先枚举0~L/2，计算$\omega_n^k$，然后枚举每一个分治块的相应位置进行蝴蝶变换。
 
IDFT
 
逆离散傅里叶变换 
 求出点值表示，再相乘之后，我们要进行插值操作。 
 通过对DFT求逆矩阵，我们直接给出以下结论： 
 将点值表示$(A_0,A_1...A_n)$转换为系数表示$(a_0,a_1...a_n)$，只需要求出： 
 
 $$a_i=\frac{\sum_{i=0}^{n-1}A_i(\omega_n^k)^{-i}}n$$  

 换句话说，我们要把DFT中的每一个

$\omega_n^k$换成

$\omega_n^{-k}$，然后除以n即可。
 
代码
 
void FFT(Complex *A,int n,int sgn)
{
    for (int i=1;i<n-1;++i)
    {
        int j=0;
        for (int t=1,k=i;t<n;t<<=1,j=((j<<1)|(k&1)),k>>=1);
        if (j>i) swap(A[i],A[j]);
    }
    for (int L=2;L<=n;L<<=1)
    {
        int L1=L>>1;
        for (int i=0;i<L1;++i)
        {
            Complex w=Complex(cos(sgn*PI*i/L1),sin(sgn*PI*i/L1));
            for (int j=i;j<n;j+=L)
            {
                int k=j+L1;
                Complex u=A[j],v=w*A[k];
                A[j]=u+v;A[k]=u-v;
            }
        }
    }
    if (sgn==-1) for (int i=0;i<n;++i) A[i]=A[i]/Complex(n,0.0);
}
 
至于是用exp还是直接两个三角函数，一个exp应该是快些，但是只能用STL的complex。如果要手写complex就只能两个三角函数。手写会快些。 
 注意：因为要对2分治，所以要强行把项数补到$2^n$，后面的项系数为0。如果是两个长度为n的多项式相乘，那么至少要开4n的空间。
 
循环卷积和线性卷积对于DFT的区别
 
线性卷积卷得到项数是2n-1，而循环卷积得到的项数是n。 
 那么我们在求点值表示的时候，线性卷积就带入2n-1个点，循环卷积就带入n个点，再相应的IDFT出来的结果就是所求的结果。
 
CZT
 
Z-变换 
 在用FFT的过程中，我们要强行把多项式补到$2^n$项。那么，FFT出来的点值表示，实际上和原来的点值表示已经不一样了。（因为n变了） 
 在求线性卷积的时候，补足$2^n$项没有什么问题。但是如果要求循环卷积，补足$2^n$位就会产生很大的问题。（定义决定的）
 
求循环卷积时，我们要保证项数不变，DFT出来的点值表示才能表示原来的循环卷积。
 
怎样保证项数不变呢？
 
考虑原来的公式 
 
 $$A_k=\sum_{j=0}^{n-1}a_j(\omega_n^k)^j$$  

 即

 $$A_k=\sum_{j=0}^{n-1}a_je^{\frac {2\pi jk}nI}$$ 
 
考虑我们把jk变换一下。 
 其实是更加复杂了 
 $jk=\frac {j^2+k^2-(j-k)^2}2$ 
 然后带入 
 
 $$A_k=\sum_{j=0}^{n-1}a_je^{\frac {\pi( j^2+k^2-(j-k)^2)}nI}$$  

 稍微变化下 

 

 $$A_k=e^{\pi\frac{k^2}nI}\sum_{j=0}^{n-1}a_je^{\frac {\pi j^2}nI}e^{-\frac{\pi (k-j)^2}nI}$$  

 设

$B_j=a_je^{\frac {\pi j^2}nI}$ 

 

$C_j=e^{-\frac{\pi j^2}nI}$ 

 那么

 $$A_k=e^{\pi\frac{k^2}nI}\sum_{j=0}^{n-1}a_jB_jC_{k-j}$$ 
 
右边是一个线性卷积，所以我们可以用一次FFT来完成一次DFT
 
注意，k-j会小于0！ 
 怎么解决这个问题呢？ 
 我们把C右移n位。那么此时 
 $C_j=e^{-\frac{\pi (j-n)^2}nI}$ 
 由于平方的存在，所以不会有问题。
 
此时C的长度变为2n。 
 卷积完之后总长3n。 
 我们需要考虑哪些部分是有用的。 
 哪些部分可以保证k-j>=0呢？ 
 即$A_{n..2n-1}$。 
 那么，最后的答案应该是：$A_k=e^{\pi\frac{k^2}nI}A_{k+n}$ 
 注意IDFT时同样要除以n。
 
这样我们就完成了项数不变的DFT。实际上，应该叫做CZT。
 
NTT
 
数论变换 
 考虑这个问题：乘起来后系数要模怎么办？
 
由于单位根是个复数，不能够模。 
 所以我们要找一个可以代替的东西。
 
考虑我们用了哪些性质，简单来说，${\omega_n^k}$是一个n阶的循环群。
 
设我们现在模的质数为$P=2^t\cdot Q+1$ 
 考虑模P意义下的原根g。
 
由于$g^{0...P-1}$两两不同，$g^P=g$，$g^{2^{t-1}Q}=-1$，所以是一个$2^t\cdot Q$阶的循环群。
 
进一步，我们可以得到$(g^Q)^{0..P-1}$是$2^t$阶的循环群。 
 让$n=2^t$即可。
 
然后我们对应着$\omega_n^k$的定义，定义$G_n^k=g^{\frac {2^tkQ}n}$。 
 显然$G^k_n$有n个，而且满足$G^k_n=G^{\frac k2}_{\frac n2}$，也满足$G_n^{k+\frac n2}=-G_n^k$
 
预处理出$g^{kQ}$代替单位根即可。
 
如果我们要模m，但是m不是$2^tQ+1$形式的质数怎么办？ 
 暴力法：在模m的情况下，每一个系数最大为 m-1，两个乘起来最多$(m-1)^2$，n个加起来最多$n(m-1)^2$。 
 我们只要凑一个超过$n(m-1)^2$的模数就好了。这样在NTT的过程中，就不会让系数发生变化。然后再去一项一项模。
 
怎么凑呢？我们用几个小的$2^tQ+1$形式的质数去做NTT，然后用中国剩余定理合并起来就好了。