首页 > 编程知识 正文

多项式的乘法法则及乘法公式,多项式乘除法则

时间:2023-05-04 19:34:03 阅读:195690 作者:2103

在上一篇文章中 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 L1 = strlen(A); int L2 = strlen(B); while(len <= 2 * L1 || len <= 2 * L2) len <<= 1; for(int i = 0; i < len; i++) { if(i < L1) a[i] = A[L1 - i - 1] - '0'; else a[i] = 0; if(i < L2) b[i] = B[L2 - i - 1] - '0'; else 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; a[k] = (u + t) % P; a[k + h / 2] = (u - t + 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] * 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) { //使用putchar()速度快很多 putchar(a[i] + '0'); flag = 0; } else if(!flag) putchar(a[i] + '0'); } 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;}


 

版权声明:该文观点仅代表作者本人。处理文章:请发送邮件至 三1五14八八95#扣扣.com 举报,一经查实,本站将立刻删除。