FFT+NTT自学笔记（大数乘积模板）

FTT学习笔记：

学习步骤

1.回顾一下复数的指数表示形式：

常规表示：
$z=a+ib；$
指数表示：
$z=re^{iθ}；$
其中，r为z的模，θ为辐角主值。
$z=r(cosθ+isinθ)=re^{iθ}$，欧拉公式

2.快速傅里叶变换要解决的问题：

计算多项式乘积（$nlogn$）
两个概念：1）离散傅里叶正变换：将多项式的系数形式转换成点值达式（复杂度（$O(nlogn$））
2）离散傅里叶反变换：将点值形式还原成多项式的系数形式（复杂度（$O(nlogn$））

3.单位复根

在一个单位圆中，将单位圆划分成n等分没个偏角度数为$θ=\frac{2*\pi}{n}$，则$e^{iθ}$即为n次单位向量，记作$W^1_n$。其余的单位向量为：$W^2_n，W^3_n，...，W^{n-1}_n$。
性质一：${(W^k_n)^n}=1$，证明：$e^{\frac{2*\pi}{n}i*kn}$=$e^{2\pi*ki}$=1
性质二：${(W^{2k}_{2n})}$=${(W^k_n)}$，证明：$e^{\frac{2*\pi}{2n}i*2k}$=$e^{\frac{2*\pi}{n}i*k}$（折半定理）
性质三：${(W^{k+\frac{n}{2}}_{n})}$=${-(W^k_n)}$，证明：$e^{\frac{2*\pi}{n}i*(k+\frac{n}{2})}$=$e^{\pi*k*i}*e^{\frac{2*\pi}{n}i*k}$=$-e^{\frac{2*\pi}{n}i*k}$=${-(W^k_n)}$（消去引理）

4.正反变换过程证明

参考大神证明

主要注意：反变换只需将$\pi$取反，重复一遍正变换过程，然后得到的多项式系数除以n即可

蝴蝶操作证明：

贴一个简洁模板：

快速傅里叶变换大数乘积简洁模板
poj2389
大数乘法模板：

#include<cstring>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<complex>
using namespace std;
#define cd complex<double>
const int maxn=2000000;//用于存放2^k个数，尽量放大
const double pi=acos(-1);
cd a[maxn],b[maxn];
int rev[maxn];

void getrev(int bit)//递推求二进制倒序后的数
{
    int len=(1<<bit);
    for(int i=0;i<len;i++)
    {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
}
void fft(cd *a,int n,int dft)
{
    //通过找规律，得到递归的最下层数的序列，然后由数组记录最下层的序列后
    //通过数组模拟往上递归过程；
    for(int i=0;i<n;i++)
    {
        if(i<rev[i])
        {
            swap(a[i],a[rev[i]]);
        }
    }
    for(int step=1;step<n;step<<=1)
    {
        cd wn=exp(cd(0,dft*pi/step));//计算本轮操作的单位复根
        for(int j=0;j<n;j+=step<<1)
        {
            cd wnk(1,0);
            for(int k=j;k<j+step;k++)
            {
                cd x=a[k],y=wnk*a[k+step];
                a[k]=x+y;//根据折半引理可得
                a[k+step]=x-y;//根据消去引理可得
                wnk*=wn;
            }
        }
    }
    if(dft==-1)//如果是逆变换
    {
        for(int i=0;i<n;i++)
        {
            a[i]/=n;
        }
    }
}
char s1[maxn],s2[maxn];
int ans[maxn];
int main()
{
    while(scanf("%s%s",s1,s2)!=EOF)
    {
        int lens1=strlen(s1);
        for(int i=0;i<lens1;i++)
        {
            a[i]=s1[lens1-i-1]-'0';
        }
        int lens2=strlen(s2);
        for(int i=0;i<lens2;i++)
        {
            b[i]=s2[lens2-i-1]-'0';
        }
        int bit=1,s=2;
        for(;(1<<bit)<lens1+lens2-1;bit++)s<<=1;
        getrev(bit);fft(a,s,1);fft(b,s,1);
        for(int i=0;i<s;i++)a[i]*=b[i];
        fft(a,s,-1);
        for(int i=0;i<s;i++)
        {
            ans[i]+=(int)(a[i].real()+0.5);
            ans[i+1]+=ans[i]/10;
            ans[i]%=10;
        }
        int i;
        for(i=lens1+lens2;!ans[i]&&i>=0;i--);
        if(i==-1)printf("0");
        for(;i>=0;i--)
        {
            printf("%d",ans[i]);
        }
        printf("\n");
    }
    return 0;
}

---

NTT学习笔记

1.了解原根

原根：如果g是p的原根，则有$g^i$ mod p ≠ $g^j$ mod p；（i≠j && 1≤i≤p-
1 && 1≤j≤p-1）
如果p是素数，设$g_n$ = $g^{\frac{p-1}{n}}$，则有一下结论：
性质一：$g_n^n=-g^{\frac{n}{2}}_n$

$$g_n^n=g^{\frac{p-1}{n}*n}_n=g^{p-1}_n=1$$
$$g^{\frac{p-1}{2}}_n=\sqrt{1}$$
由原根的定义可知当 $p-1$ ≠ $\frac{p-1}{2}$， $g^{\frac{p-1}{2}}_n$ ≠ 1，则 $g^{\frac{p-1}{2}}_n$=-1;性质一得证
性质二： $g_n^{2k}=g^k_{\frac{n}{2}}$
$$g_n^{2k}=g^{\frac{p-1}{n}*2k}_n$$
$$g^{k}_{\frac{n}{2}}=g^{\frac{p-1}{n}*2k}_{\frac{n}{2}}$$
得证；（以上所有结论都要模p）

以上两个性质正好和fft的单位复根作为插值是完全一样的，因此可以用原根模数作为插值，正反变换过程不变，但是未用到浮点数，因此精度提高。

2.关于n的选值问题

在fft中，取大于等于n的2的幂作物最终的n
在ntt中，通常选$a*2^k+1$的素数
通常选则

$$p=1004535809=479*{2^{21}}+1$$
因为它的最小正原根是3
UOJ 模数 $998244353=2^{23}×7×17+1$，最小正原根也是 3
快速数论变换大数乘积简洁模板

poj2389

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define ll long long
using namespace std;
const int mod=479*(1<<21)+1;
const int maxn=2000000;
int rev[maxn];
ll a[maxn],b[maxn];
ll qm(ll a,ll b)
{
    if(b<0)
    {
        b=-b;
        a=qm(a,mod-2);
    }
    a=a%mod;
    ll ans=1;
    while(b)
    {
        if(b&1)
        {
            ans=ans*a%mod;
        }
        b>>=1;
        a=a*a%mod;
    }
    return ans;
}

void getrev(int bit)
{
    for(int i=0; i<(1<<bit); i++)
    {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
}
void ntt(ll *a,int n,int dft)
{
    for(int i=0; i<n; i++)
    {
        if(i<rev[i])
        {
            swap(a[i],a[rev[i]]);
        }
    }
    for(int step=1; step<n; step<<=1)
    {
        ll wn;
        wn=qm(3,dft*(mod-1)/(step*2));//即3^((p-1)/n)
        for(int j=0; j<n; j+=step<<1)
        {
            ll wnk=1;
            for(int k=j; k<j+step; k++)
            {
                //通通注意mod
                ll x=a[k]%mod;
                ll y=(wnk*a[step+k])%mod;
                a[k]=(x+y)%mod;
                a[k+step]=((x-y)%mod+mod)%mod;
                wnk=(wnk*wn)%mod;
            }
        }
    }
    if(dft==-1)
    {
        int ni=qm(n,mod-2);
        for(int i=0; i<n; i++)
        {
            a[i]=a[i]*ni%mod;//注意mod
        }
    }
}
char s1[maxn],s2[maxn];
int ans[maxn];
int main()
{
    scanf("%s%s",s1,s2);
    int lens1=strlen(s1),lens2=strlen(s2);
    int bit=1,s=2;
    for(; (1<<bit)<(lens1+lens2-1); bit++)
        s<<=1;

    getrev(bit);

    for(int i=0; i<lens1; i++)
    {
        a[i]=s1[lens1-i-1]-'0';
    }
    for(int j=0; j<lens2; j++)
    {
        b[j]=s2[lens2-j-1]-'0';
    }
    ntt(a,s,1);
    ntt(b,s,1);
    for(int i=0; i<s; i++)
        a[i]=(a[i]*b[i])%mod;
    ntt(a,s,-1);
    for(int i=0; i<s; i++)
    {
        ans[i]+=a[i];
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
    }
    int i;
    for(i=s; !ans[i]&&i>=0; i--);
    if(i==-1)
        printf("0\n");
    for(; i>=0; i--)
    {
        printf("%d",ans[i]);
    }
    printf("\n");

    return 0;
}
/*123423
34564356546*/