题解 P1919 【【模板】A*B Problem升级版（FFT快速傅里叶）】

· · 题解

鉴于毒瘤鰰加强了数据,在此发一篇关于加强数据的NTT题解

一、如何优化高精乘

众所周知,朴素的高精乘是 \Theta(n^2) 的,对于超过 10^{10000} 的数据,要尝试去优化ta
一个比较简单的方法是压位,就是用一个int或者long long连续存储多位,以此优化常数,然而...


于是 a,b\le 10^{65536} => a,b\le 10^{1000000} ,就必须优化复杂度来做这道题了。怎么优化呢?

首先,我们可以通过位值原理将一个数转化为n次函数(多项式),
如\quad3579=3x^3+5x^2+7x+9\quad,其中 x=10。

我们称这种我们一般用来表示n此函数的方法称为系数表示法。
不难发现,该多项式的乘法就是不带进位的高精乘 (加法同理)

此外,用小学知识可知,用 (n+1) 个点可唯一确定一个 n 次函数,
我们定义这种表示方法为点值表示法。
于是设第一个多项式A可转化为 (x_0,a_0),(x_1,a_1)\cdots(x_n,a_n)，
第二个多项式B可转化为\quad\quad\quad(x_0,b_0),(x_1,b_1)\cdots(x_n,b_n)，
我们就可以用 \Theta(n) 的复杂度把它们乘起来,于是答案就是

(x_0,a_0\times b_0),(x_1,a_1\times b_1)\cdots(x_n,a_n\times b_n)

Q:可是代入这些值也要 \Theta(n^2) 呀,而且显而易见会爆long long, 这种做法不应更劣吗?
A:不急,我们接下来就讲。

二、傅里叶的优化方法

鉴于这一情况,大佬傅里叶找到了一组特殊的值在特定意义下代入,以此来缩短值域,优化常数和精度

我们把 将系数表示法 通过代入他找到的值 转化为点值表示法 称为DFT,
再把 由DFT转化成的点值表示法 重新转化为系数表示法 称为IDFT

如何让代入的值达到这种效果呢? 我们希望代入的值满足以下条件：
1.代入的值互不相同(显而易见,如果相同就相当于同一个点)
2.代入的值的乘方在一个范围内转圈圈。(这个比较重要,对值域和精度都有限制)

于是我们就有两种做法来限制:取模(NTT)和复数运算(FFT)。其中复数运算涉及涉及高中内容,
相关定理对于初中萌新较为难理解,而且实数运算有精度限制,所以这里我们用取模来限制,即NTT

三、NTT的前置芝士

1.乘法逆元(鰰叕在模板上面卡了常)
P3811 【模板】乘法逆元
P5431 【模板】乘法逆元2
2.原根
学过乘法逆元的朋友应该都知道费马小定理: \text{若}p\text{为质数,则}a^{p-1}\equiv 1\quad(\text{mod }p)
那么对于一个质数 p ,如果有一个数 a ,使得 a^0,a^1,\cdots a^{p-2} 在\mod p意义下互不相等，
,我们就称a为\mod m的原根,通常用字母 g 表示

换个方式理解,原根的乘方相当于把1\sim p-1全部遍历了一遍 然后回到了原点,刚好满足条件1和条件2

四、NTT加速多项式乘法

1.单位根
为了方便分治,我们先将多项式补足到2^k,
然后将形如 q\times 2^k+1 的质数作为模数,让计算在其模意义下进行

NTT常用的模数有998244353=119\times 2^{23}+1,1004535809=479\times 2^{21},二者的原根都为3
如果需要更多的模数,可以看这里

然后我们代入的数就是

g^0,g^q,g^{2q},g^{3q}\cdots g^{(2^k-1)q}(\text{mod }q\times 2^k+1)

于是在模意义下代入它们就得到了DFT,
而IDFT就是就是按顺序把代入得到的值作为系数,依次代入这些数的逆元,
最后乘上2^k的逆元就可以转回来了

为方便称呼,我们设2^k=n,w^x_n=g^{xq},此时w^x_n我们称它为单位根

于是就可以得到一些性质：

$\quad$2.$w^{2x}_{2n}\equiv w^x_n(\text{mod } q\times 2^k+1)$(相当于从$2^k$中提取一个2到q) $\quad$3.若$x<\frac{n}{2}$,则$w^x_n\equiv -w^{x+\frac{n}{2}}_n(\text{mod }q\times 2^k+1) $\qquad\qquad 1\equiv -w^{\frac{n}{2}}_n(\text{mod }q\times 2^k+1) \qquad\qquad\because (w^{\frac{n}{2}}_n)^2=w^n_n=1,$且由原根的性质$,w^{\frac{n}{2}}_n\ne w^0_n=1 透过这些性质，就可以用分治来做了 设要代入的多项式为$F(x)=a_0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1}

把下标为偶数的放一边,下标为奇数的放一边,得到

F(x)=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots a_{n-1}x^{n-2})

于是设A(x)=a_0+a_2x^1+\cdot+a_{n-2}x^{\frac{n}{2}-1},B(x)=a_1+a_3x^1+\cdot+a_{n-1}x^{\frac{n}{2}-1}
所以F(x)=A(x^2)+x\times B(x^2)
代入w^x_n,得到F(w^x_n)=A((w^x_n)^2)+w^x_n\times B((w^x_n)^2)=A(w^x_{\frac{n}{2}})+w^x_n\times B(w^x_{\frac{n}{2}})
这时,若我们代入w^{x+\frac{n}{2}}_n,就可以得到

F(w^{x+\frac{n}{2}}_n)=A((w^{x+\frac{n}{2}}_n)^2)+w^{x+\frac{n}{2}}_n\times B((w^{x+\frac{n}{2}}_n)^2)=A(w^x_{\frac{n}{2}})-w^x_n\times B(w^x_{\frac{n}{2}})

这时,我们发现,我们在扫过F(w^0_n\sim w^{\frac{n}{2}-1}_n)时,可以顺带得到F(w^{\frac{n}{2}}_n\sim w^{n-1}_n),也就可以\Theta(n\log n)分治了

五、FFT/NTT的常数优化

1.非递归写法
我们在划分A(x),B(x)的时候,是按下标的奇偶性来分的,也就是看下标二进制的最后一位
于是整体递归下来,就是先看下标二进制的最后一位,再看倒数第二位,再...
于是最终递归到的位置就是下标二进制翻转过来,可以预处理好

2.蝴蝶变换：就地进行FFT\NTT
如果已经顺序存储了A(1\sim w^{\frac{n}{4}}_{\frac{n}{2}}),B(1\sim w^{\frac{n}{4}}_{\frac{n}{2}})，A在B的前面
对于要迭代到的w^x_n和w^{x+\frac{n}{2}}_n，可以同时进行迭代,
将w^x_n放到原A(w^x_{\frac{n}{2}})的位置,将w^{x+\frac{n}{2}}放到原B(w^x_{\frac{n}{2}})的位置,保证顺序不变

六、最终代码

#include<stdio.h>
#include<string.h>
typedef unsigned char byte;
typedef unsigned int word;
typedef unsigned long long ull;
typedef long long ll;
const word mod=1004535809,size=21;
//NTT模数,结果不超过(1<<size) 
char io[(1<<size)+1];//输入输出 
word a[1<<size],b[1<<size];//多项式/高精 
word realid[1<<size],root[1<<size],inverse[1<<size];
//迭代后系数所在的位置,单位根及其逆元 
word i,id,floor;
ll num1,num2;
char *top;
//循环变量 
inline ll pow(ll a,ll b){//快速幂 
    register ll ans=1;
    for(;b;b>>=1){
        if(b&1) (ans*=a)%=mod;
        (a*=a)%=mod;
    }
    return ans;
}
inline void loading(){//预处理迭代后位置,单位根及其逆元 
    root[0]=inverse[0]=1;
    num1=pow(3,mod>>size);
    num2=pow(num1,mod-2);
    for(i=1;i<1<<size;i++){
        root[i]=num1*root[i-1]%mod;
        inverse[i]=num2*inverse[i-1]%mod;
        for(id=i,floor=0;floor<size;floor++,id>>=1)
            realid[i]=realid[i]<<1|(id&1);
    }
}
inline void read(){
//直接用fread,然后指针前移读入数据,效率较高,写法较简单 
//因fread特性,DEV-CPP下不能以数字结尾(不停读入最后一个字符),但linux下可以 
    top=io+(1<<size);
    fread(io+1,1,1<<size,stdin);
    while('0'>*top||*top>'9') top--;
    for(i=0;'0'<=*top&&*top<='9';top--,i++)
        a[realid[i]]=*top-'0';//直接放到迭代后的位置 
    while('0'>*top||*top>'9') top--;    
    for(i=0;'0'<=*top&&*top<='9';top--,i++)
        b[realid[i]]=*top-'0';//b同理 
} 
inline void DFT(){ //非迭代版DFT 
    for(floor=0;floor<size;floor++)
        for(i=0;i<1<<size;i+=(1<<(floor+1)))
            for(id=0;id<1<<floor;id++){
                num1=a[i+id];//蝴蝶变换 
                num2=a[i+id+(1<<floor)];
                (num2*=root[id<<size-floor-1])%=mod;
                a[i+id]=(num1+num2)%mod;//放回原位 
                a[i+id+(1<<floor)]=(num1+mod-num2)%mod;

                num1=b[i+id];//b同理 
                num2=b[i+id+(1<<floor)];
                (num2*=root[id<<size-floor-1])%=mod;
                b[i+id]=(num1+num2)%mod;
                b[i+id+(1<<floor)]=(num1+mod-num2)%mod;
            }
}
inline void IDFT(){
    for(floor=0;floor<21;floor++)//与DFT相同 
        for(i=0;i<0x200000;i+=1<<floor+1)
            for(id=0;id<1<<floor;id++){
                num1=root[i+id];
                num2=root[i+id+(1<<floor)];
                (num2*=inverse[id<<20-floor])%=mod;//乘上单位根逆元 
                root[i+id]=(num1+num2)%mod;
                root[i+id+(1<<floor)]=(num1+mod-num2)%mod;
            }
}
inline void write(){//fwrite输出,同输入
    num2=pow(1<<size,mod-2);
    top=io+(1<<size);
    num1=0;
    for(i=0;i<1<<size;i++){
        num1+=num2*root[i]%mod;//最后要乘上n的逆元 
        top--;
        *top=num1%10+'0';
        num1/=10;
    }
    while(*top=='0'&&top!=io+(1<<size)) top++;
    if(top==io+(1<<size)) putchar('0');
    else fwrite(top,1,io+(1<<size)-top,stdout);
}
int main(){
    loading();
    read(); 
    DFT();
    for(i=0;i<1<<size;i++)
        root[realid[i]]=(ll)(a[i])*b[i]%mod;//点值相乘 
    IDFT();
    write(); 
    return 0;
}

END、数据范围

设NTT模数为q\times 2^k+1,那么多项式长度不能超过2^k

设原本相乘多项式的系数不超过n,长度为P,则相乘后的系数不超过n^2P

如果要使得到的结果在模意义下仍旧是唯一的,除了扩大NTT模数,也可以取多个模数做NTT,
最后将得到的结果用中国剩余定理合并(这也就是任意模数NTT的方法)