快速数论变换（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）实现。

 

  1 #include <iostream>  
  2 #include <string.h>  
  3 #include <stdio.h>  
  4   
  5 using namespace std;  
  6 typedef long long LL;  
  7   
  8 const int N = 1 << 18;  
  9 const int P = (479 << 21) + 1;  
 10 const int G = 3;  
 11 const int NUM = 20;  
 12   
 13 LL  wn[NUM];  
 14 LL  a[N], b[N];  
 15 char A[N], B[N];  
 16   
 17 LL quick_mod(LL a, LL b, LL m)  
 18 {  
 19     LL ans = 1;  
 20     a %= m;  
 21     while(b)  
 22     {  
 23         if(b & 1)  
 24         {  
 25             ans = ans * a % m;  
 26             b--;  
 27         }  
 28         b >>= 1;  
 29         a = a * a % m;  
 30     }  
 31     return ans;  
 32 }  
 33   
 34 void GetWn()  
 35 {  
 36     for(int i=0; i<NUM; i++)  
 37     {  
 38         int t = 1 << i;  
 39         wn[i] = quick_mod(G, (P - 1) / t, P);  
 40     }  
 41 }  
 42   
 43 void Prepare(char A[], char B[], LL a[], LL b[], int &len)  
 44 {  
 45     len = 1;  
 46     int len_A = strlen(A);  
 47     int len_B = strlen(B);  
 48     while(len <= 2 * len_A || len <= 2 * len_B) len <<= 1;  
 49     for(int i=0; i<len_A; i++)  
 50         A[len - 1 - i] = A[len_A - 1 - i];  
 51     for(int i=0; i<len - len_A; i++)  
 52         A[i] = '0';  
 53     for(int i=0; i<len_B; i++)  
 54         B[len - 1 - i] = B[len_B - 1 - i];  
 55     for(int i=0; i<len - len_B; i++)  
 56         B[i] = '0';  
 57     for(int i=0; i<len; i++)  
 58         a[len - 1 - i] = A[i] - '0';  
 59     for(int i=0; i<len; i++)  
 60         b[len - 1 - i] = B[i] - '0';  
 61 }  
 62   
 63 void Rader(LL a[], int len)  
 64 {  
 65     int j = len >> 1;  
 66     for(int i=1; i<len-1; i++)  
 67     {  
 68         if(i < j) swap(a[i], a[j]);  
 69         int k = len >> 1;  
 70         while(j >= k)  
 71         {  
 72             j -= k;  
 73             k >>= 1;  
 74         }  
 75         if(j < k) j += k;  
 76     }  
 77 }  
 78   
 79 void NTT(LL a[], int len, int on)  
 80 {  
 81     Rader(a, len);  
 82     int id = 0;  
 83     for(int h = 2; h <= len; h <<= 1)  
 84     {  
 85         id++;  
 86         for(int j = 0; j < len; j += h)  
 87         {  
 88             LL w = 1;  
 89             for(int k = j; k < j + h / 2; k++)  
 90             {  
 91                 LL u = a[k] % P;  
 92                 LL t = w * (a[k + h / 2] % P) % P;  
 93                 a[k] = (u + t) % P;  
 94                 a[k + h / 2] = ((u - t) % P + P) % P;  
 95                 w = w * wn[id] % P;  
 96             }  
 97         }  
 98     }  
 99     if(on == -1)  
100     {  
101         for(int i = 1; i < len / 2; i++)  
102             swap(a[i], a[len - i]);  
103         LL Inv = quick_mod(len, P - 2, P);  
104         for(int i = 0; i < len; i++)  
105             a[i] = a[i] % P * Inv % P;  
106     }  
107 }  
108   
109 void Conv(LL a[], LL b[], int n)  
110 {  
111     NTT(a, n, 1);  
112     NTT(b, n, 1);  
113     for(int i = 0; i < n; i++)  
114         a[i] = a[i] * b[i] % P;  
115     NTT(a, n, -1);  
116 }  
117   
118 void Transfer(LL a[], int n)  
119 {  
120     int t = 0;  
121     for(int i = 0; i < n; i++)  
122     {  
123         a[i] += t;  
124         if(a[i] > 9)  
125         {  
126             t = a[i] / 10;  
127             a[i] %= 10;  
128         }  
129         else t = 0;  
130     }  
131 }  
132   
133 void Print(LL a[], int n)  
134 {  
135     bool flag = 1;  
136     for(int i = n - 1; i >= 0; i--)  
137     {  
138         if(a[i] != 0 && flag)  
139         {  
140             printf("%d", a[i]);  
141             flag = 0;  
142         }  
143         else if(!flag)  
144             printf("%d", a[i]);  
145     }  
146     puts("");  
147 }  
148   
149 int main()  
150 {  
151     GetWn();  
152     while(scanf("%s%s", A, B)!=EOF)  
153     {  
154         int len;  
155         Prepare(A, B, a, b, len);  
156         Conv(a, b, len);  
157         Transfer(a, len);  
158         Print(a, len);  
159     }  
160     return 0;  
161 }

 

posted @ 2017-08-05 00:54  HuaZhang  阅读(6519)  评论(0)  编辑  收藏  举报