码迷,mamicode.com
首页 > 其他好文 > 详细

多项式乘法运算终极版

时间:2014-09-05 18:15:26      阅读:229      评论:0      收藏:0      [点我收藏+]

标签:style   blog   http   color   os   io   ar   strong   for   

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

换来求多项式的乘法。可以发现它是利用了单位复根的特殊性质,大大减少了运算,但是这种做法是对复数系数的矩阵

加以处理,每个复数系数的实部和虚部是一个正弦及余弦函数,因此大部分系数都是浮点数,我们必须做复数及浮点数

的计算,计算量会比较大,而且浮点数的计算可能会导致误差增大。

 

今天,我将来介绍另一种计算多项式乘法的算法,叫做快速数论变换(NTT),在离散正交变换的理论中,已经证明在

复数域内,具有循环卷积特性的唯一变换是DFT,所以在复数域中不存在具有循环卷积性质的更简单的离散正交变换。

因此提出了以数论为基础的具有循环卷积性质的快速数论变换

 

回忆复数向量,其离散傅里叶变换公式如下

 

   bubuko.com,布布扣

 

离散傅里叶逆变换公式为

 

   bubuko.com,布布扣

 

今天的快速数论变换(NTT)是在bubuko.com,布布扣上进行的,在快速傅里叶变换(FFT)中,通过bubuko.com,布布扣次单位复根来运算的,即满

bubuko.com,布布扣bubuko.com,布布扣,而对于快速数论变换来说,则是可以将bubuko.com,布布扣看成是bubuko.com,布布扣的等价,这里bubuko.com,布布扣是模素数bubuko.com,布布扣

的原根(由于bubuko.com,布布扣是素数,那么原根一定存在)。即

 

        bubuko.com,布布扣

 

所以综上,我们得到数论变换的公式如下

 

    bubuko.com,布布扣

 

数论变换的逆变换公式为

 

    bubuko.com,布布扣

 

这样就把复数对应到一个整数,之后一切都是在bubuko.com,布布扣系统内考虑。

 

上述数论变换(NTT)公式中,要求bubuko.com,布布扣是素数且bubuko.com,布布扣必须是bubuko.com,布布扣的因子。由于bubuko.com,布布扣经常是2的方幂,所以可以构造形

bubuko.com,布布扣的素数。通常来说可以选择bubuko.com,布布扣费马素数,这样的变换叫做费马数数论变换

 

这里我们选择bubuko.com,布布扣bubuko.com,布布扣,这样得到模bubuko.com,布布扣的原根值为bubuko.com,布布扣

 

 

题目: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 LL N = 1 << 18;
const LL P = (479 << 21) + 1;
const LL 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_A = strlen(A);
    int len_B = strlen(B);
    for(int i=0; i<len_A; i++)
        A[N - 1 - i] = A[len_A - 1 - i];
    for(int i=0; i<N - len_A; i++)
        A[i] = '0';
    for(int i=0; i<len_B; i++)
        B[N - 1 - i] = B[len_B - 1 - i];
    for(int i=0; i<N - len_B; i++)
        B[i] = '0';
    for(int i=0; i<N; i++)
        a[N - 1 - i] = A[i] - '0';
    for(int i=0; i<N; i++)
        b[N - 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 % P * (a[k + h / 2] % P) % P;
                a[k] = (u % P + t % P) % P;
                a[k + h / 2] = ((u % P - t % P) % P + P) % P;
                w = w % P * 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)
    {
        Prepare(A, B, a, b);
        Conv(a, b, N);
        Transfer(a, N);
        Print(a, N);
    }
    return 0;
}

 

多项式乘法运算终极版

标签:style   blog   http   color   os   io   ar   strong   for   

原文地址:http://blog.csdn.net/acdreamers/article/details/39026505

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!