NC212291. FFT快速傅里叶
描述
输入描述
第一行一个正整数n。 第二行描述一个位数为n的正整数x。 第三行描述一个位数为n的正整数y。
输出描述
输出一行，即x*y的结果。
示例1
输入:
1 3 4
输出:
12
说明:
n<=60000C++(g++ 7.5.0) 解法, 执行用时: 39ms, 内存消耗: 7360K, 提交时间: 2022-08-14 20:26:51
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
const double PI = acos(-1.0);
struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x;
y = _y;
}
Complex operator-(const Complex &b) const {
return Complex(x - b.x, y - b.y);
}
Complex operator+(const Complex &b) const {
return Complex(x + b.x, y + b.y);
}
Complex operator*(const Complex &b) const {
return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
/*
* 进行 FFT 和 IFFT 前的反置变换
* 位置 i 和 i 的二进制反转后的位置互换
*len 必须为 2 的幂
*/
void change(Complex y[], int len) {
int i, j, k;
for (int i = 1, j = len / 2; i < len - 1; i++) {
if (i < j) std::swap(y[i], y[j]);
// 交换互为小标反转的元素，i<j 保证交换一次
// i 做正常的 + 1，j 做反转类型的 + 1，始终保持 i 和 j 是反转的
k = len / 2;
while (j >= k) {
j = j - k;
k = k / 2;
}
if (j < k) j += k;
}
}
/*
* 做 FFT
*len 必须是 2^k 形式
*on == 1 时是 DFT，on == -1 时是 IDFT
*/
void fft(Complex y[], int len, int on) {
change(y, len);
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
for (int j = 0; j < len; j += h) {
Complex w(1, 0);
for (int k = j; k < j + h / 2; k++) {
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn;
}
}
}
if (on == -1) {
for (int i = 0; i < len; i++) {
y[i].x /= len;
}
}
}
const int MAXN = 200020;
Complex x1[MAXN], x2[MAXN];
char str1[MAXN / 2], str2[MAXN / 2];
int sum[MAXN];
int main() {
int n;std::cin>>n;
scanf("%s%s", str1, str2);
int len1 = strlen(str1);
int len2 = strlen(str2);
int len=1;
while (len < len1 * 2 || len < len2 * 2) len <<= 1;
for (int i = 0; i < len1; i++) x1[i] = Complex(str1[len1 - 1 - i] - '0', 0);
for (int i = len1; i < len; i++) x1[i] = Complex(0, 0);
for (int i = 0; i < len2; i++) x2[i] = Complex(str2[len2 - 1 - i] - '0', 0);
for (int i = len2; i < len; i++) x2[i] = Complex(0, 0);
fft(x1, len, 1);
fft(x2, len, 1);
for (int i = 0; i < len; i++) x1[i] = x1[i] * x2[i];
fft(x1, len, -1);
for (int i = 0; i < len; i++) sum[i] = int(x1[i].x + 0.5);
for (int i = 0; i < len; i++) {
sum[i + 1] += sum[i] / 10;
sum[i] %= 10;
}
len = len1 + len2 - 1;
while (sum[len] == 0 && len > 0) len--;
for (int i = len; i >= 0; i--) printf("%c", sum[i] + '0');
printf("\n");
return 0;
}
C++14(g++5.4) 解法, 执行用时: 55ms, 内存消耗: 5880K, 提交时间: 2020-09-23 19:27:33
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=4e5+5;
const double PI=acos(-1.0);
int n,m;
struct cp{
double x,y;
inline cp operator +(const cp b)const{return (cp){x+b.x,y+b.y};}
inline cp operator -(const cp b)const{return (cp){x-b.x,y-b.y};}
inline cp operator *(const cp b)const{return (cp){x*b.x-y*b.y,x*b.y+y*b.x};}
}a[maxn],b[maxn];
int lenth=1,Reverse[maxn],ans[maxn];
void init(int x){
int tim=0;lenth=1;
while(lenth<=x)lenth<<=1,tim++;
for(int i=0;i<=lenth;i++){
a[i]={0,0},b[i]={0,0},Reverse[i]=0;
}
for(int i=0;i<lenth;i++)Reverse[i]=(Reverse[i>>1]>>1)|((i&1)<<(tim-1));
}
void FFT(cp *A,const int fla){
for(int i=0;i<lenth;i++)if(i<Reverse[i])swap(A[i],A[Reverse[i]]);
for(int i=1;i<lenth;i<<=1){
cp w=(cp){cos(PI/i),fla*sin(PI/i)};
for(int j=0;j<lenth;j+=(i<<1)){
cp K=(cp){1,0};
for(int k=0;k<i;k++,K=K*w){
cp x=A[j+k],y=A[j+k+i]*K;
A[j+k]=x+y;
A[j+k+i]=x-y;
}
}
}
}
int main(){
int p;cin>>p;
string x,y;
while(cin>>x>>y){
n=x.size()-1,m=y.size()-1;
reverse(x.begin(),x.end());
reverse(y.begin(),y.end());
init(n+m);
for(int i=0;i<=n;i++)a[i].x=x[i]-'0';
for(int i=0;i<=m;i++)b[i].x=y[i]-'0';
FFT(a,1),FFT(b,1);
for(int i=0;i<lenth;i++)a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<=n+m;i++)ans[i]=(int)(a[i].x/lenth+0.5);
int k=0;
string s;
for(int i=0;i<=n+m;i++){
s+=(char)((ans[i]+k)%10+'0');
k=(ans[i]+k)/10;
}
while(k>0)s+=k%10+'0',k/=10;
reverse(s.begin(),s.end());
k=0;
while(s[k]=='0'&&k+1<s.size())k++;
while(k<s.size())cout<<s[k++];
cout<<'\n';
}
}
C++ 解法, 执行用时: 234ms, 内存消耗: 189564K, 提交时间: 2022-07-13 15:05:30
//#include <bits/stdc++.h>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<complex>
using namespace std;
typedef long long ll;
typedef complex<double> comp;
const int MAXN=1e6+5;
const comp I(0,1);
const double PI=acos(-1);
comp A[MAXN*3],B[MAXN*3],tmp[MAXN*3],ans[MAXN*3];//数组开大一些
int rev[MAXN*3];
void fft(comp F[],int N,int sgn=1){
int bit=__lg(N);
for(int i=0;i<N;++i){// 蝴蝶变换
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
if(i<rev[i]){
swap(F[i],F[rev[i]]);
}
}
for(int l=1;l<N;l<<=1){//枚举合并前的区间长度
comp step=exp(sgn*PI/l*I);
for(int i=0;i<N;i+=l*2){//遍历起始点
comp cur(1,0);
for(int k=i;k<i+l;++k){
comp g=F[k],h=F[k+l]*cur;
F[k]=g+h,F[k+l]=g-h;
cur*=step;
}
}
}
}
inline int read(){
int ans = 0;
char c = getchar();
while (!isdigit(c))
c = getchar();
while (isdigit(c))
{
ans = ans * 10 + c - '0';
c = getchar();
}
return ans;
}
int o[MAXN];
int main()
{
int n;
cin>>n;
string s1,s2;
cin>>s1>>s2;
for(int i=0;i<n;i++){
A[i]=s1[i]-'0';
B[i]=s2[i]-'0';
}
int N = 1 << __lg(n + n + 1) + 1;
fft(A, N), fft(B, N);
for (int i=0;i<N;++i){
ans[i]=A[i]*B[i];
}
fft(ans,N,-1);
for(int i=0;i<n+n-1;++i){
o[i]=int(ans[i].real()/N+0.1);
}
// for(int i=0;i<n+n-1;++i){
// cout<<o[i]<<" ";
// }
// cout<<endl;
int t=0;
for(int i=n+n-2;i>=0;i--){
o[i]+=t;
t=o[i]/10;
o[i]%=10;
}
if(t!=0){
cout<<t;
}
for(int i=0;i<n+n-1;++i){
cout<<o[i];
}
return 0;
}
/*
5 5
1 1 1 1 1
1 1 1 1 1
*/
Python3 解法, 执行用时: 343ms, 内存消耗: 8180K, 提交时间: 2022-09-22 21:31:03
n=int(input()) x=int(input()) y=int(input()) print(x*y)
pypy3 解法, 执行用时: 208ms, 内存消耗: 25788K, 提交时间: 2023-07-21 20:36:38
n=int(input()) a=int(input()) c=int(input()) print(a*c)