您的位置:首页 > 其它

【UOJ#34】 多项式乘法(FFT && NTT)

2017-01-30 20:59 525 查看
记录一个菜逼的成长。。

题目链接

这是一道模板题。。

FFT版本:

#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long LL;
const double PI = acos(-1);
const int INF = 0x3f3f3f3f;
struct complex
{
double r,i;
complex(double _r = 0.0,double _i = 0.0)
{
r = _r; i = _i;
}
complex operator +(const complex &b)
{
return complex(r+b.r,i+b.i);
}
complex operator -(const complex &b)
{
return complex(r-b.r,i-b.i);
}
complex operator *(const complex &b)
{
return complex(r*b.r-i*b.i,r*b.i+i*b.r);
}
};
///雷德算法--倒位序
void Rader(complex F[],int len)
{
for( int i = 1,j = len >> 1; i < len-1; i++ ){
if(i < j)swap(F[i],F[j]);
int k = len >> 1;
while(j >= k){
j -= k;
k >>= 1;
}
if(j < k)j += k;
}
}
/**
* 做FFT
* len必须为2^k形式,
* on==1时是DFT,on==-1时是IDFT
**/
void fft(complex F[],int len,int dft)
{
Rader(F,len);
for( int h = 2; h <= len; h <<= 1 ){
complex wn(cos(-dft*2*PI/h),sin(-dft*2*PI/h));
for( int j = 0; j < len; j += h ){
complex w(1,0);
for( int k = j; k < j + (h >> 1); k++ ){
complex u = F[k];
complex t = w*F[k+(h>>1)];
F[k] = u + t;
F[k+(h>>1)] = u - t;
w = w * wn;
}
}
}
if(dft == -1)
for( int i = 0; i < len; i++ )
F[i].r /= len;
}
void Conv(complex a[],complex b[],int len)
{
fft(a,len,1);
fft(b,len,1);
for( int i = 0; i < len; i++ )
a[i] = a[i] * b[i];
fft(a,len,-1);
}
const int maxn = 500000 + 10;
int n,m;
int aa[maxn],bb[maxn];
complex a[maxn],b[maxn];
int main()
{
while(~scanf("%d%d",&n,&m)){
for( int i = 0; i <= n; i++ )
scanf("%d",aa+i);
for( int i = 0; i <= m; i++ )
scanf("%d",bb+i);
int len = 1;
while(len < (n+1)*2 || len < (m+1)*2)len <<= 1;
for( int i = 0; i <= n; i++ )
a[i] = complex(aa[i],0);
for( int i = n+1; i < len; i++ )
a[i] = complex(0,0);
for( int i = 0; i <= m; i++ )
b[i] = complex(bb[i],0);
for( int i = m+1; i < len; i++ )
b[i] = complex(0,0);
Conv(a,b,len);
for(int i = 0; i < n + m + 1; i++ ){
if(!i)printf("%d",(int)(a[i].r+0.5));
else printf(" %d",(int)(a[i].r+0.5));
}
puts("");
}
return 0;
}


NTT版本:

可以参考

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

using namespace std;
typedef long long LL;

const int maxn = 1 << 18;
const int P = (479 << 21) + 1;
const int G = 3;
const int NUM = 20;
LL  wn[NUM];
LL  a[m
e6ab
axn], b[maxn];
char A[maxn], B[maxn];
LL PowerMod(LL a, LL b, LL c)
{
LL ans = 1;
a = a % c;
while(b>0)
{
if(b % 2 == 1)
ans = (ans * a) % c;
b = b/2;
a = (a * a) % c;
}
return ans;
}
void GetWn()
{
for(int i = 0; i < NUM; i++)
{
int t = 1 << i;
wn[i] = PowerMod(G, (P - 1) / t, P);
}
}
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 = PowerMod(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);
}
int main()
{
GetWn();
int n,m;
while(~scanf("%d%d",&n,&m)){
for( int i = 0; i <= n; i++ )scanf("%d",a+i);
for( int i = 0; i <= m; i++ )scanf("%d",b+i);
int len = 1;
while(len < 2*(n+1) || len < 2*(m+1))len <<= 1;
for( int i = n + 1; i < len; i++ )a[i] = 0;
for( int i = m + 1; i < len; i++ )b[i] = 0;
Conv(a,b,len);
for( int i = 0; i < n + m + 1; i++ ){
if(!i)printf("%d",a[i]);
else printf(" %d",a[i]);
}
puts("");
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签:  fft 多项式乘法 ntt