初識FFT和NTT

看了很久文檔，覺得自己隻學會了套模闆的能力，了解的代碼是怎麼寫，還有一點原理，看完現在來推一下原理估計又不會了！

初識FFT和NTT

學這個的原因是因為codechef的一道題目，可惜現在還是沒有解決

初識FFT和NTT

，誰會了求教點選打開連結

看了 ACdream 的部落格，覺得不夠詳細，而且看了之後根本看不懂代碼裡面寫的什麼鬼，之後找了一份題解

然後在百度文庫找到了一篇講得很詳細的文檔，看玩總算了解了那麼一點點！

初識FFT和NTT

ps：這裡我隻能寫一點代碼要做什麼，原理不敢亵渎！

1、首先是反轉置換，将一個2的n次幂長度的分成兩個，奇數是一組，偶數是一組，然後一次遞歸下去，直到長度為1，這樣每個元素的位置會去到他二進制反轉之後變成是十進制的位置，做置換可以采用rader算法線性實作！一下引用上面連結部落格的一張圖：

初識FFT和NTT

其實為什麼這樣想想也很明顯，将位置pos編号換成二進制，初始化ans_pos=0，ans依次做位運算右移，低位為0的時候走向二叉樹的左兒子，否則是右兒子，而ans_pos依次左移，向二叉樹左兒子走的時候末尾補0，否則補1，這樣ans_pos的二進制就和剛開始ans的二進制是相反的

rader算法代碼如下：

///雷德算法，2^M=len，将第i位的數與“i的二進制反轉之後的位”的數交換
void rader(complex *F,int len)
{
    int j=len/2;///模拟二進制反轉進位的的位置
    for(int i=1;i<len-1;i++)
    {
        if(i<j)swap(F[i],F[j]);///該出手時就出手
        int k=len/2;
        while(j>=k)
        {
            j-=k;
            k>>=1;
        }
        if(j<k)j+=k;
    }
}

然後抄了份模闆，具體解釋就看上面連結ppt的解釋，非常詳細：

const int N = 500005;
const double pi = acos(-1.0);
struct complex{///複數類
    double real,imag;
    complex (double r=0.0,double i=0.0)
    {
        real=r;
        imag=i;
    }
    complex operator + (const complex &x)
    {
        return complex(x.real+real,x.imag+imag);
    }
    complex operator - (const complex &x)
    {
        return complex(real-x.real,imag-x.imag);
    }
    complex operator * (const complex &x)
    {
        return complex(x.real*real-x.imag*imag,imag*x.real+real*x.imag);
    }
}vara[N],varb[N];
void FFT(complex *F,int len,int t)
{
    rader(F,len);
    for(int h=2;h<=len;h<<=1)///L from 1 to M,2^M=len
    {
        ///J_product_two_power_of_M_subtract_L
        ///M_subtract_L_product_h_equals_to_len
        complex wn(cos(t*2*pi/h),sin(t*2*pi/h));///公比
        for(int j=0;j<len;j+=h)
        {
            complex E(1,0);///螺旋因子
            for(int k=j;k<j+h/2;k++)
            {
                complex u=F[k];///蝶型操作
                complex v=E*F[k+h/2];
                F[k]=u+v;///前半部分
                F[k+h/2]=u-v;///後半部分
                E=E*wn;
            }
        }
    }
    if(t==-1)///IDFT
    for(int i=0;i<len;i++)
        F[i].real/=len;
}

NTT也是搞一發模闆：

ACdream 點選打開連結

點選打開連結這個一份好的解釋看完可以了解一二了

大數運算NTT：

const int N = 400005;
const int g=3;
const long long MOD=(479<<21)+1;
int len;
char A[N],B[N];
long long a[N],b[N],qp[30];
long long q_pow(long long x,long long y,long long P)
{
    long long ans=1;
    while(y>0)
    {
        if(y&1)ans=ans*x%P;
        x=x*x%P;
        y>>=1;
    }
    return ans;
}
void init()
{
    len=1;
    for(int i=0;i<21;i++)
    {
        int t=1<<i;
        qp[i]=q_pow(g,(MOD-1)/t,MOD);
    }
    int len1=strlen(A);
    int len2=strlen(B);
    while(len<=2*len1||len<=2*len2)len<<=1;
    for(int i=0;i<len1;i++)a[len1-i-1]=A[i]-'0';
    for(int i=len1;i<len;i++)a[i]=0;
    for(int i=0;i<len2;i++)b[len2-i-1]=B[i]-'0';
    for(int i=len2;i<len;i++)b[i]=0;
}
///雷德算法，2^M=len，将第i位的數與“i的二進制反轉之後的位”的數交換
void rader(long long F[],int len)
{
    int j=len/2;///模拟二進制反轉進位的的位置
    for(int i=1;i<len-1;i++)
    {
        if(i<j)swap(F[i],F[j]);///該出手時就出手
        int k=len/2;
        while(j>=k)
        {
            j-=k;
            k>>=1;
        }
        if(j<k)j+=k;
    }
}
void NTT(long long F[],int len,int t)
{
    int id=0;
    rader(F,len);
    for(int h=2;h<=len;h<<=1)
    {
        id++;
        for(int j=0;j<len;j+=h)
        {
            long long E=1;///原根次幂
            for(int k=j;k<j+h/2;k++)
            {
                long long u=F[k];///蝶型操作
                long long v=(E*F[k+h/2])%MOD;
                F[k]=(u+v)%MOD;///前半部分
                F[k+h/2]=((u-v)%MOD+MOD)%MOD;///後半部分
                E=(E*qp[id])%MOD;
            }
        }
    }
    //p(F);
    if(t==-1)///插值
    {
        for(int i=1;i<len/2;i++)swap(F[i],F[len-i]);///i+lne-i=i;
        long long inv=q_pow(len,MOD-2,MOD);///逆元
        for(int i=0;i<len;i++)F[i]=(F[i]%MOD*inv)%MOD;
    }
    //p(F);
}
void work()///卷積，點乘，插值
{
    NTT(a,len,1);
    NTT(b,len,1);
    for(int i=0;i<len;i++)
        a[i]=(a[i]*b[i])%MOD;
    NTT(a,len,-1);
}
void pushup()///進位
{
    for(int i=0;i<len;i++)
    {
        if(a[i]>=10)
        {
            a[i+1]+=a[i]/10;
            a[i]=a[i]%10;
        }
    }
}
void print()///輸出
{
    int high=0;
    for(int i=len-1;i>=0;i--)
    {
        if(a[i])
        {
            high=i;
            break;
        }
    }
    for(int i=high;i>=0;i--)printf("%lld",a[i]);
    puts("");
}

hdu1402

點選打開連結

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

using namespace std;

const int N = 400005;
const double pi = acos(-1.0);
const int g=3;
//const int len=1<<18;
const long long MOD=(479<<21)+1;
int len;
char A[N],B[N];
long long a[N],b[N],qp[30];
void p(long long a[])
{
    for(int i=0;i<len;i++)printf("%lld ",a[i]);puts("");
}
long long q_pow(long long x,long long y,long long P)
{
    long long ans=1;
    while(y>0)
    {
        if(y&1)ans=ans*x%P;
        x=x*x%P;
        y>>=1;
    }
    return ans;
}
void init()
{
    len=1;
    for(int i=0;i<21;i++)
    {
        int t=1<<i;
        qp[i]=q_pow(g,(MOD-1)/t,MOD);
    }
    int len1=strlen(A);
    int len2=strlen(B);
    while(len<2*len1||len<2*len2)len<<=1;
    for(int i=0;i<len1;i++)a[len1-i-1]=A[i]-'0';
    for(int i=len1;i<len;i++)a[i]=0;
    for(int i=0;i<len2;i++)b[len2-i-1]=B[i]-'0';
    for(int i=len2;i<len;i++)b[i]=0;
}
void rader(long long F[],int len)
{
    int j=len/2;
    for(int i=1;i<len-1;i++)
    {
        if(i<j)swap(F[i],F[j]);
        int k=len/2;
        while(j>=k)
        {
            j-=k;
            k>>=1;
        }
        if(j<k)j+=k;
    }
}
void NTT(long long F[],int len,int t)
{
    int id=0;
    rader(F,len);
    for(int h=2;h<=len;h<<=1)
    {
        id++;
        for(int j=0;j<len;j+=h)
        {
            long long E=1;
            for(int k=j;k<j+h/2;k++)
            {
                long long u=F[k];
                long long v=(E*F[k+h/2])%MOD;
                F[k]=(u+v)%MOD;
                F[k+h/2]=((u-v)%MOD+MOD)%MOD;
                E=(E*qp[id])%MOD;
            }
        }
    }
    //p(F);
    if(t==-1)
    {
        for(int i=1;i<len/2;i++)swap(F[i],F[len-i]);
        long long inv=q_pow(len,MOD-2,MOD);
        for(int i=0;i<len;i++)F[i]=(F[i]%MOD*inv)%MOD;
    }
    //p(F);
}
void work()
{
    NTT(a,len,1);
    NTT(b,len,1);
    for(int i=0;i<len;i++)
        a[i]=(a[i]*b[i])%MOD;
    NTT(a,len,-1);
}
void pushup()
{
    for(int i=0;i<len;i++)
    {
        if(a[i]>=10)
        {
            a[i+1]+=a[i]/10;
            a[i]=a[i]%10;
        }
    }
}
void print()
{
    int high=0;
    for(int i=len-1;i>=0;i--)
    {
        if(a[i])
        {
            high=i;
            break;
        }
    }
    for(int i=high;i>=0;i--)printf("%lld",a[i]);
    puts("");
}
int main()
{
    while(~scanf("%s%s",A,B))
    {
        init();
        work();
        pushup();
        print();
    }
    return 0;
}

初識FFT和NTT

繼續閱讀

【數學】HDU 1222 Wolf and Rabbit

多屬性決策模型詳解（matlab）

2020年高教社模組化國賽真題B題--穿越沙漠

網際網路面試題之趣味數學題

多元函數微分學強化班

【數學】MIT線性代數總結

ZOJ 1007 Numerical Summation of a Series

BWT (Burrows–Wheeler_transform) 解碼分析

Day07-數學基礎-經濟問題(DataWhale)

HDOJ 3664 Permutation Counting / UVALive 5092 DP

STFT filter bankSTFT filter bank

有關拉格朗日反演的擴充形式的證明

相位相關算法

【CF 125D】 Two progressions 劃分等差數列

Codeforces Round #319 (Div. 2) B C D

2017多校訓練Contest4: 1003 Counting Divisors hdu6069