快速數論變換（NTT）

轉自ACdreamers （http://blog.csdn.net/acdreamers/article/details/39026505）

 

 

在上一篇文章中 http://blog.csdn.net/acdreamers/article/details/39005227 介紹了用快速傅里葉變

換來求多項式的乘法。可以發現它是利用了單位復根的特殊性質，大大減少了運算，但是這種做法是對復數系數的矩陣

加以處理，每個復數系數的實部和虛部是一個正弦及余弦函數，因此大部分系數都是浮點數，我們必須做復數及浮點數

的計算，計算量會比較大，而且浮點數的計算可能會導致誤差增大。

 

今天，我將來介紹另一種計算多項式乘法的算法，叫做快速數論變換（NTT），在離散正交變換的理論中，已經證明在

復數域內，具有循環卷積特性的唯一變換是DFT，所以在復數域中不存在具有循環卷積性質的更簡單的離散正交變換。

因此提出了以數論為基礎的具有循環卷積性質的快速數論變換。

 

回憶復數向量，其離散傅里葉變換公式如下

 

   

 

離散傅里葉逆變換公式為

 

   

 

今天的快速數論變換（NTT）是在上進行的，在快速傅里葉變換（FFT）中，通過次單位復根來運算的，即滿

足的，而對於快速數論變換來說，則是可以將看成是的等價，這里是模素數

的原根（由於是素數，那么原根一定存在）。即

 

        

 

所以綜上，我們得到數論變換的公式如下

 

    

 

而數論變換的逆變換公式為

 

    

 

這樣就把復數對應到一個整數，之后一切都是在系統內考慮。

 

上述數論變換（NTT）公式中，要求是素數且必須是的因子。由於經常是2的方冪，所以可以構造形

如的素數。通常來說可以選擇為費馬素數，這樣的變換叫做費馬數數論變換。

 

這里我們選擇，，這樣得到模的原根值為。

 

 

題目：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1028

 

分析：題目意思就是大數相乘，此處用快速數論變換（NTT）實現。

 

#include <iostream>  
#include <string.h>  
#include <stdio.h>  
  
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 len_A = strlen(A);  
    int len_B = strlen(B);  
    while(len <= 2 * len_A || len <= 2 * len_B) len <<= 1;  
    for(int i=0; i<len_A; i++)  
        A[len - 1 - i] = A[len_A - 1 - i];  
    for(int i=0; i<len - len_A; i++)  
        A[i] = '0';  
    for(int i=0; i<len_B; i++)  
        B[len - 1 - i] = B[len_B - 1 - i];  
    for(int i=0; i<len - len_B; i++)  
        B[i] = '0';  
    for(int i=0; i<len; i++)  
        a[len - 1 - i] = A[i] - '0';  
    for(int i=0; i<len; i++)  
        b[len - 1 - i] = 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) % P;  
                a[k] = (u + t) % P;  
                a[k + h / 2] = ((u - t) % P + 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] % P * 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)  
        {  
            printf("%d", a[i]);  
            flag = 0;  
        }  
        else if(!flag)  
            printf("%d", a[i]);  
    }  
    puts("");  
}  
  
int main()  
{  
    GetWn();  
    while(scanf("%s%s", A, B)!=EOF)  
    {  
        int len;  
        Prepare(A, B, a, b, len);  
        Conv(a, b, len);  
        Transfer(a, len);  
        Print(a, len);  
    }  
    return 0;  
}