Dylan 收录于 OI

2018-03-13 2023-04-05 约 900 字预计阅读 2 分钟 - 次阅读 - 条评论

目录

快速傅里叶变换

多项式乘法

例题：luogu3803

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
/*
luogu3803
FFT
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
using namespace std;
template<typename T>
void read(T& w)
{
	char rev;int f=1;
	for(rev=getchar();(rev<48||rev>57)&&rev!='-';rev=getchar());
	if(rev=='-')rev=getchar(),f=-1;
	for(w=0;rev>=48&&rev<=57;rev=getchar())w=w*10+rev-48;
	w*=f;
}
template<typename T>
void write(T w,char c=0)
{
	if(w<0)putchar('-'),w=-w;
	if(w>9)write(w/10);
	putchar(w%10+48);
	if(c)putchar(c);
}
template<typename T>
void Swap(T& a,T& b){T t=a;a=b;b=t;}
const double pi=3.14159265358979323846264338;
template<typename T>
struct COMPLEX
{
	T a,b;
	COMPLEX(){a=b=0;}
	COMPLEX(T aa,T bb):a(aa),b(bb){}
	T& real(){return a;}
	T& imag(){return b;}
};

template<typename T>
COMPLEX<T> operator+(const COMPLEX<T>& a,const COMPLEX<T>& b){return COMPLEX<T>(a.a+b.a,a.b+b.b);}
template<typename T>
COMPLEX<T> operator-(const COMPLEX<T>& a,const COMPLEX<T>& b){return COMPLEX<T>(a.a-b.a,a.b-b.b);}
template<typename T>
COMPLEX<T> operator*(const COMPLEX<T>& a,const COMPLEX<T>& b){return COMPLEX<T>(a.a*b.a-a.b*b.b,a.a*b.b+a.b*b.a);}
template<typename T>
COMPLEX<T> operator*=(COMPLEX<T>& a,const COMPLEX<T>& b){a=COMPLEX<T>(a.a*b.a-a.b*b.b,a.a*b.b+a.b*b.a);}
typedef COMPLEX<double> C;
const int maxn=2*int(1e6+3e5);
int n,m,l,lg,rev[maxn];
C a[maxn],b[maxn],ans[maxn];

void FFT(C* A,int f=1)
{
	for(int i=0;i<=lg-1;++i)
		if(i<rev[i])Swap(A[i],A[rev[i]]);
	for(int i=1;i<=lg-1;i<<=1)
	{
		C wn(cos(pi/i),f*sin(pi/i));
		for(int p=i<<1,j=0;j<=lg-1;j+=p)
		{
			C w(1,0);
			for(int k=0;k<=i-1;++k,w*=wn)
			{
				C u=A[j+k],v=w*A[j+k+i];
				A[j+k]=u+v;
				A[j+k+i]=u-v;
			}
		}
	}
}
#define IFFT(A) FFT(A,-1)
void input()
{
	read(n),read(m);
	for(int i=0;i<=n;++i)read(a[i].real());
	for(int i=0;i<=m;++i)read(b[i].real());
}
void initRev()
{
	for(lg=1;lg<=n+m;lg<<=1)++l;//l=log2(n+m)+1,n=2^l
	for(int i=0;i<=lg-1;++i)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));//¶þ½øÖÆÄæÐò×ª»» 
}
void solve()
{
	FFT(a),FFT(b);
	for(int i=0;i<=lg;++i)
		ans[i]=a[i]*b[i];
	IFFT(ans);
	for(int i=0;i<=n+m;++i)
		ans[i].real()/=lg;
}
void output()
{
	for(int i=0;i<=n+m;++i)
		write(int(ans[i].real()+0.5),i==n+m?'\n':' ');
}
int main()
{
	freopen("luogu3803_1.in","r",stdin);
	//freopen("luogu3803.out","w",stdout);
	input();
	initRev();
	solve();
	output();
	return 0;
}

快速高精乘

例题：luogu1919

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
/*
luogu1919
FFT
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
using namespace std;
int getint()
{
	char r;
	do{r=getchar();}
	while(r<48||r>57);
	return r-48;
}
template<typename T>
void write(T w,char c=0)
{
	if(w<0)putchar('-'),w=-w;
	if(w>9)write(w/10);
	putchar(w%10+48);
	if(c)putchar(c);
}
template<typename T>
void Swap(T& a,T& b){T t=a;a=b;b=t;}
const double pi=3.14159265358979323846264338;
const double eps=1e-9;
template<typename T>
struct COMPLEX
{
	T a,b;
	COMPLEX(){a=b=0;}
	COMPLEX(T aa,T bb):a(aa),b(bb){}
	T& real(){return a;}
	T& imag(){return b;}
};

template<typename T>
COMPLEX<T> operator+(const COMPLEX<T>& a,const COMPLEX<T>& b){return COMPLEX<T>(a.a+b.a,a.b+b.b);}
template<typename T>
COMPLEX<T> operator-(const COMPLEX<T>& a,const COMPLEX<T>& b){return COMPLEX<T>(a.a-b.a,a.b-b.b);}
template<typename T>
COMPLEX<T> operator*(const COMPLEX<T>& a,const COMPLEX<T>& b){return COMPLEX<T>(a.a*b.a-a.b*b.b,a.a*b.b+a.b*b.a);}
template<typename T>
COMPLEX<T> operator*=(COMPLEX<T>& a,const COMPLEX<T>& b){return a=a*b;}
template<typename T>
COMPLEX<T> operator+=(COMPLEX<T>& a,const COMPLEX<T>& b){return a=a+b;}
typedef COMPLEX<double> C;
const int maxn=4*60003;
int n,m,l,lg,rev[maxn];
C a[maxn],b[maxn],c[maxn];
int ans[maxn];
void FFT(C* A,int f=1)
{
	for(int i=0;i<=lg-1;++i)
		if(i<rev[i])Swap(A[i],A[rev[i]]);
	for(int i=1;i<=lg-1;i<<=1)
	{
		C wn(cos(pi/i),f*sin(pi/i));
		for(int p=i<<1,j=0;j<=lg-1;j+=p)
		{
			C w(1,0);
			for(int k=0;k<=i-1;++k,w*=wn)
			{
				C u=A[j+k],v=w*A[j+k+i];
				A[j+k]=u+v;
				A[j+k+i]=u-v;
			}
		}
	}
}
#define IFFT(A) FFT(A,-1)
void input()
{
	scanf("%d",&n);
	--n;
	for(int i=n;i>=0;--i)a[i].real()=getint();
	for(int i=n;i>=0;--i)b[i].real()=getint();
}
void initRev()
{
	for(lg=1;lg<=2*n;lg<<=1)++l;//l=log2(2n)+1,n=2^l
	for(int i=0;i<=lg-1;++i)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));//¶þ½øÖÆÄæÐò×ª»» 
}
void solve()
{
	FFT(a),FFT(b);
	for(int i=0;i<=lg;++i)
		c[i]=a[i]*b[i];
	IFFT(c);
	for(int i=0;i<=2*n;++i)
		c[i].real()/=lg;
	for(int i=0;i<=2*n;++i)
	{
		ans[i]+=int(c[i].real()+0.5);
		if(ans[i]>9)ans[i+1]+=ans[i]/10,ans[i]%=10;
	}
}
void output()
{
	for(int i=2*n+1,flag=0;i>=0;--i)
	{
		int x=ans[i];
		if(x&&!flag)flag=1;
		if(flag)write(x,i==0?'\n':0);
	}
}
int main()
{
	freopen("luogu1919_t.in","r",stdin);
	//freopen("luogu1919.out","w",stdout);
	input();
	initRev();
	solve();
	output();
	return 0;
}