DFT/FFT/NTT

在Seal库和HElib库中都用到了NTT技术，用于加快多项式计算，而NTT又是FFT的优化，FFT又来自于DFT，现在具体学习一下这三个技术！

基础概念#

名词区分#

1、DFT：离散傅立叶变换
2、FFT：快速傅立叶变换
3、NTT：快速数论变换
4、MTT：NTT的扩展
5、多项式卷积：多项式乘法
6、根据多项式的系数表示法求点值表示法的过程叫做“求值”；根据点值表示法求系数表示法的过程称为“插值”
7、求一个多项式的乘法，即求卷积，先通过傅立叶变换对系数表示法的多项式进行求值运算，其复杂度 $O (n l o g^{n})$ ，然后在 $O (n)$ 的时间内点值相乘，在进行插值运算。

8、如果选取单位复根作为求值点，则可以对系数向量进行离散傅立叶变换（DFT），得到相应的点值表示；同样可以通过对点值进行逆DFT运算，获得相应的系数向量。DFT和逆DFT时间复杂度均为 $O (n l o g^{n})$ 。

复数#

定义#

我们知道，一个复数可以这样表示： $a + b i$ ，a和b是实数，其中 $i$ 叫做虚数单位，复数域是目前已知最大的域。
在复平面中，x轴代表实数，y轴（除原点外的点）代表虚数，从原点(0,0)到(a,b)的向量表示复数 $a + b i$

模长：从原点(0,0)到(a,b)的距离，即 $\sqrt{a^{2} + b^{2}}$
幅角：假设以逆时针为正方向，从x轴正半轴到已知向量的转角的有向角叫做幅角

运算#

1、加法
在复数平面，复数可以表示为向量，因为复数的加法和向量的加法相同。
2、乘法
几何定义：复数相乘密，模长相乘，幅角相加
代数定义：

(a + b i) * (c + d i) = a c + a d i + b c i + b d i^{2} + a c + a d i + b c i - b d = (a c - b d) + (b c + a d) i

单位根#

在复数平面上，以原点为圆心，1为半径做圆，所得的圆叫做单位圆，以圆点为起点，圆的n等分为终点，做n个向量，设幅角为正且最小的向量对应的复数为w_n$，称为n次单位根。

根据复数乘法的运算法则，其余n-1个复数为 $w_{n}^{2}, . . ., w_{n}^{n}$ ，注意 $w_{n}^{0} = w_{n}^{n} = 1$ （对应复平面上以x轴为正方向的向量）

如何计算呢？
由欧拉公式解决 $w_{n}^{k} = c o s (k * 2 π / n) + i s i n (k * 2 π / n)$
例如：向量AB表示的是复数为4次单位根

n次单位根的幅角为周角的 $1 / n$ 。

在代数中，若 $z_{n} = 1$ ，我们把z称为n次单位根。具体请参考：n次单位根（n-th unit root）

单位根的性质#

1、 $w_{n}^{k} = c o s (k * 2 π / n) + i s i n (k * 2 π / n)$

2、【相消引理】 $w_{d n}^{d k} = w_{n}^{k}$
证明：以d=2为例

3、【折半引理】 $w_{n}^{k + n / 2} = - w_{n}^{k}$
证明：

4、 $w_{n}^{0} = w_{n}^{k} = 1$

5、 $w_{n}^{n - i} = 共轭 (w_{n}^{i})$

6、 $w_{n}^{n + i} = w_{n}^{i}$

多项式系数表示法#

设 $A (x)$ 表示一个d次多项式，则 $A (x) = a_{1} + a_{2} * x + . . ., + a_{d} * x^{d}$
利用这种方法计算多项式卷积复杂度为 $O (d^{2})$ ，其实就是直接对应相乘（暴力）。
例如： $A (x) = 1 + 2 x + x^{2}$ , $B (x) = 1 - 2 x + x^{2}$

A (x) B (x) = (+ 2 x + x^{2}) (1 - 2 x + x^{2}) = 1 - 2 x^{2} + x^{4}

多项式点值表示法#

将n个值x带入多项式，会得到d各不同的值y，则该多项式被这n个点值 $(x_{1}, y_{1}), . . ., (x_{d}, y_{x})$ 唯一确定，其中 $\sum_{j = 1}^{d} a_{j} * x_{i}^{j}$
而利用点值法计算多项式卷积复杂度也为 $O (d^{2})$ 。（选点 $O (d)$ ，每次计算 $O (d)$ ）
例如上面的多项式用点值法表示： $A (x) = [(- 2, 1), (- 1, 0), (0, 1), (1, 4), (2, 9)], B (x) = [(- 2, 9), (- 1, 4), (0, 1), (1, 0), (2, 1)]$ ，则

C (x) = [(- 2, 9), (- 1, 0), (0, 1), (1, 0), (2, 9)]

即有这个5个点就可以唯一确定一个4次多项式，而两两相乘的复杂度为 $O (d)$

引理1： $(d + 1)$ 个点值可以唯一确定一个d 阶多项式

因此，我们可以将一个系数多项式转换为一个点值多项式，然而进行复杂度为 $O (d)$ 的乘法，再将结果的点值多项式恢复回系数多项式。

但是：如果我们采用下面这种矩阵形式计算点值的话【选点】，那么由系数转为点值的复杂度也为 $O (d^{2})$ 。

接下来考虑对其优化：
1、对于系数表示法，每个点的系数都固定，优化困难
2、对于点值表示法，可以用FFT来解决！

DFT#

已知 $A (x)$ 的系数为 $(a_{0}, a_{1}, . . ., a_{n - 1})$ ，对于 $k = 0, 1, . . ., n - 1$ ，定义：

y_{k} = A (w_{n}^{k}) = \sum_{i = 0}^{n - 1} a_{i} w_{n}^{k i}

其中向量 $y = (y_{0}, y_{1}, . . ., y_{n - 1})$ 是系数向量 $a = (a_{0}, a_{1}, . . ., a_{n - 1})$ 的离散傅立叶变换，记 $y = D F T_{n} (a)$ ，复杂度为 $O (n^{2})$
而使用下面的FFT方法，可以在 $O (n l o g^{n})$ 时间内求出 $D F T_{n} (a)$

FFT#

用于加速系数多项式到点值多项式的运算！
首先观察下面多项式：

例如： $F (x) = x^{2}$ ，有对称性 $F (- x) = F (x)$ ，相当于确定了一个点相当于确定两个点。
同理又如 $F (x^{3})$ ，有性质 $F (- x) = - F (x)$ ，也是确定了一个点相当于确定了两个点。

所以对于有奇偶行的多项式，只需要找到原本一半的点就可以得到这个多现实了。
基于以上想法，假如有下面多项式：

把 $P_{e}$ 和 $P_{o}$ 分别看作两个多项式，也就是对于一个点 $x_{i}$ ，我们只要计算出 $P_{e} (x_{i}^{2})$ 和 $P_{o} (x_{i}^{2})$ ，就可以得到 $P (x_{i})$ 和 $P (- x_{i})$ ，而且 $P_{e}$ 和 $P_{o}$ 还可以进一步拆分为奇偶两部分！

假设原本我们需要n个点 $【 \pm x_{1}, \pm x_{1} 2, . . ., \pm x_{n / 2} 】$ 就能确定一个 $n - 1$ 阶的多项式。现在变成了求 $P_{e} (x)$ 和 $P_{o} (x)$ 在 $x_{1}^{2}, x_{2}^{2}, . . ., x_{n / 2}^{2}$ 上面的点值【n/2个点】。
那如果这n/2个点两两之间满足 $x_{i}^{2} = - x_{j}^{2}$ ，则就可以进一步拆分为一半了，就可以将原本的复杂度 $O (d^{2})$ 降为 $O (d l o g^{d})$ 。这里可以看出FFT用到了分治思想。

问题是，$x_1^2,x_22,...,x_{n/2}^2并不满足两两互为相反数。由此使用n次单位根，选用n个n次单文根

$[w^{0}, w^{1}, . . ., w^{n - 1}]$ 。

这样，两个点平方后依旧互为相反数！

可以看出，将以一个n个点的求值问题转换为求n/2个点，在转换为求n/4个点，以此迭代，从而达到 $O (d l o g^{d})$ 。将上述思想转换为为代码如下：

FFT的逆#

如何从点值多项式变为系数多项式呢？

对于点值计算，

实际上就是一个矩阵的乘法：

将点换为n个n次单位根，则矩阵变为：

其中中间的范德蒙德矩阵就成了一个DFT矩阵。
有了正向（系数到点值）的矩阵变换，求逆向（点值到系数）就是对上面矩阵求逆即可：

即：

从上面可以看出，FFT是将 $w$ 作为点值传入，IFFT就是将 $1 / n * w^{- 1}$ 作为点值传入：

程序#

下面程序用FFT计算两个大数乘
题目：http://acm.hdu.edu.cn/showproblem.php?pid=1402

#include <iostream>
#include <string.h>
#include <stdio.h>
#include <math.h>

using namespace std;
const int N = 500005;
const double PI = acos(-1.0);

struct Virt
{
    double r, i;

    Virt(double r = 0.0,double i = 0.0)
    {
        this->r = r;
        this->i = i;
    }

    Virt operator + (const Virt &x)
    {
        return Virt(r + x.r, i + x.i);
    }

    Virt operator - (const Virt &x)
    {
        return Virt(r - x.r, i - x.i);
    }

    Virt operator * (const Virt &x)
    {
        return Virt(r * x.r - i * x.i, i * x.r + r * x.i);
    }
};

//雷德算法--倒位序
void Rader(Virt F[], int len)
{
    int j = len >> 1;
    for(int i=1; i<len-1; i++)
    {
        if(i < j) swap(F[i], F[j]);
        int k = len >> 1;
        while(j >= k)
        {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}

//FFT实现
void FFT(Virt F[], int len, int on)
{
    Rader(F, len);
    for(int h=2; h<=len; h<<=1) //分治后计算长度为h的DFT
    {
        Virt wn(cos(-on*2*PI/h), sin(-on*2*PI/h));  //单位复根e^(2*PI/m)用欧拉公式展开
        for(int j=0; j<len; j+=h)
        {
            Virt w(1,0);            //旋转因子
            for(int k=j; k<j+h/2; k++)
            {
                Virt u = F[k];
                Virt t = w * F[k + h / 2];
                F[k] = u + t;     //蝴蝶合并操作
                F[k + h / 2] = u - t;
                w = w * wn;      //更新旋转因子
            }
        }
    }
    if(on == -1)
        for(int i=0; i<len; i++)
            F[i].r /= len;
}

//求卷积
void Conv(Virt a[],Virt b[],int len)
{
    FFT(a,len,1);
    FFT(b,len,1);
    for(int i=0; i<len; i++)
        a[i] = a[i]*b[i];
    FFT(a,len,-1);
}

char str1[N],str2[N];
Virt va[N],vb[N];
int result[N];
int len;

void Init(char str1[],char str2[])
{
    int len1 = strlen(str1);
    int len2 = strlen(str2);
    len = 1;
    while(len < 2*len1 || len < 2*len2) len <<= 1;

    int i;
    for(i=0; i<len1; i++)
    {
        va[i].r = str1[len1-i-1] - '0';
        va[i].i = 0.0;
    }
    while(i < len)
    {
        va[i].r = va[i].i = 0.0;
        i++;
    }
    for(i=0; i<len2; i++)
    {
        vb[i].r = str2[len2-i-1] - '0';
        vb[i].i = 0.0;
    }
    while(i < len)
    {
        vb[i].r = vb[i].i = 0.0;
        i++;
    }
}

void Work()
{
    Conv(va,vb,len);
    for(int i=0; i<len; i++)
        result[i] = va[i].r+0.5;
}

void Export()
{
    for(int i=0; i<len; i++)
    {
        result[i+1] += result[i]/10;
        result[i] %= 10;
    }
    int high = 0;
    for(int i=len-1; i>=0; i--)
    {
        if(result[i])
        {
            high = i;
            break;
        }
    }
    for(int i=high; i>=0; i--)
        printf("%d",result[i]);
    puts("");
}

int main()
{
    while(~scanf("%s%s",str1,str2))
    {
        Init(str1,str2);
        Work();
        Export();
    }
    return 0;
}

NTT#

在FFT中，我们需要用到复数，复数虽然很神奇，但是它也有自己的局限性——需要用double类型计算，精度太低，那有没有什么东西能够代替复数且解决精度问题呢？
这个东西，叫原根

阶#

若a,p互素，且p>1，对于 $a^{n} m o d p = 1$ 满足最小的n，叫做a模p的阶，记 $δ_{p} (a)$ .
例如：

δ_{7} (2) = 3

其中：
$2^{1} m o d 7 = 2$
$2^{2} m o d 7 = 4$
$2^{3} m o d 7 = 1$

原根#

设p是正整数，a是整数，若 $δ_{p} (a)$ 等于 $ϕ (p)$ ，则a为模p的一个原根。

例如：
$δ_{7} (3) = 6 = ϕ (7)$ ，所以3是模7的一个原根。

原根的个数不唯一

1、若模数p有原根，那么它一定有 $ϕ (ϕ (p)) 个原根$
2、若p为素数，原根一定存在，假设g是P的一个原根，那么 $g^{i} m o d p (1 < g < p, 0 < i < p)$ 的结果两两不同
简单的说，就是

g^{i} m o d p \neq g^{j} m o d p, (1 < i \neq j < p - 1)

3、那如何求一个质数的原根呢？
对于指数p， $p_{i}$ 是p-1的因子，若 $g^{p - 1 / p_{i}} （ m o d p ）$ 恒成立，则g是p的原根。

下面就是为什么原根可以代替单位根计算？
因为原根具有和单位根相同的性质，FFT中，用到了单位根的四条性质，原根也满足这四条性质：

最终可以得到:

w_{n} = g^{p - 1 / n} m o d p

然后只需将FFT中的 $w_{n}$ 替换掉，就是NTT。即：

综上，NTT的变换为：

这里P是素数且N必须是P-1的因子；由于N是2的方幂，所以可构造 $P = c {.2}^{k} + 1$ 的素数。
通常p取998244353，它的原根为3。

程序#

使用NTT，计算两个大数乘

#include <iostream>
#include <string.h>
#include <stdio.h>
#include <ctime>
using namespace std;
typedef long long LL;

const int N = 1 << 18;
const int P = (479 << 21) + 1;
const int G = 3;
const int NUM = 20;

LL  wn[NUM];
LL  a[N], b[N];
char A[N], B[N];

LL quick_mod(LL a, LL b, LL m)
{
    LL ans = 1;
    a %= m;
    while(b)
    {
        if(b & 1)
        {
            ans = ans * a % m;
            b--;
        }
        b >>= 1;
        a = a * a % m;
    }
    return ans;
}

void GetWn()
{
    for(int i = 0; i < NUM; i++)
    {
        int t = 1 << i;
        wn[i] = quick_mod(G, (P - 1) / t, P);
    }
}

void Prepare(char A[], char B[], LL a[], LL b[], int &len)
{
    len = 1;
    int L1 = strlen(A);
    int L2 = strlen(B);
    while(len <= 2 * L1 || len <= 2 * L2) len <<= 1;
    for(int i = 0; i < len; i++)
    {
        if(i < L1) a[i] = A[L1 - i - 1] - '0';
        else a[i] = 0;
        if(i < L2) b[i] = B[L2 - i - 1] - '0';
        else b[i] = 0;
    }
}

void Rader(LL a[], int len)
{
    int j = len >> 1;
    for(int i = 1; i < len - 1; i++)
    {
        if(i < j) swap(a[i], a[j]);
        int k = len >> 1;
        while(j >= k)
        {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}

void NTT(LL a[], int len, int on)
{
    Rader(a, len);
    int id = 0;
    for(int h = 2; h <= len; h <<= 1)
    {
        id++;
        for(int j = 0; j < len; j += h)
        {
            LL w = 1;
            for(int k = j; k < j + h / 2; k++)
            {
                LL u = a[k] % P;
                LL t = w * a[k + h / 2] % P;
                a[k] = (u + t) % P;
                a[k + h / 2] = (u - t + P) % P;
                w = w * wn[id] % P;
            }
        }
    }
    if(on == -1)
    {
        for(int i = 1; i < len / 2; i++)
            swap(a[i], a[len - i]);
        LL inv = quick_mod(len, P - 2, P);
        for(int i = 0; i < len; i++)
            a[i] = a[i] * inv % P;
    }
}

void Conv(LL a[], LL b[], int n)
{
    NTT(a, n, 1);
    NTT(b, n, 1);
    for(int i = 0; i < n; i++)
        a[i] = a[i] * b[i] % P;
    NTT(a, n, -1);
}

void Transfer(LL a[], int n)
{
    int t = 0;
    for(int i = 0; i < n; i++)
    {
        a[i] += t;
        if(a[i] > 9)
        {
            t = a[i] / 10;
            a[i] %= 10;
        }
        else t = 0;
    }
}

void Print(LL a[], int n)
{
    bool flag = 1;
    for(int i = n - 1; i >= 0; i--)
    {
        if(a[i] != 0 && flag)
        {
            //使用putchar()速度快很多
            putchar(a[i] + '0');
            flag = 0;
        }
        else if(!flag)
            putchar(a[i] + '0');
    }
    puts("");
}


int main()
{
    GetWn();
//71992652622957199265262295622895175513333762898450963102326778447440803476281886006800109449932463374706102636321647845385139107625719926526229571992652622956228951755133337628984509631023267784474408034762818860068001094499324633747061026363216478453851391076257199265262295719926526229562289517551333376289845096310232677844744080347628188600680010944993246337470610263632164784538513910762571992652622957199265262295622895175513333762898450963102326778447440803476281886006800109449932463374706102636321647845385139107625719926526229571992652622956228951755133337628984509631023267784474408034762818860068001094499324633747061026363216478453851391076257199265262295719926526229562289517551333376289845096310232677844744080347628188600680010944993246337470610263632164784538513910762571992652622957199265262295622895175513333762898450963102326778447440803476281886006800109449932463374706102636321647845385139107625
    while(scanf("%s %s", A, B) != EOF)
    {
        int len;
        clock_t start_time = clock();//计时开始
        Prepare(A, B, a, b, len);
        Conv(a, b, len);
        Transfer(a, len);
        cout << "elapsed time:" << 1000*double(clock() - start_time) / CLOCKS_PER_SEC
             << 'ms' << endl;
        Print(a, len);
    }
    return 0;
}