G42 快速傅里叶变换 FFT算法 高精度乘法

视频链接：https://www.bilibili.com/video/BV1EA411D7yD/

Luogu P1919 【模板】A*B Problem 升级版

// 递归版 6.6s
#include <iostream>
#include <cstring>
#include <algorithm>
#include<cmath>
using namespace std;
const int N=3e6;
const double PI=acos(-1);

struct complex{
  double x, y;
  complex operator+(const complex& t)const{
    return {x+t.x, y+t.y};}
  complex operator-(const complex& t)const{
    return {x-t.x, y-t.y};}
  complex operator*(const complex& t)const{
    return {x*t.x-y*t.y, x*t.y+y*t.x};}
}A[N], B[N];
char s1[N], s2[N];
int ans[N];

void FFT(complex A[],int n,int op){
  if(n==1)return;
  complex A1[n/2],A2[n/2];
  for(int i=0;i<n/2;++i)
    A1[i]=A[i*2], A2[i]=A[i*2+1];
  FFT(A1,n/2,op); FFT(A2,n/2,op);
  complex w1({cos(2*PI/n),sin(2*PI/n)*op});
  complex wk({1,0});
  for(int i=0;i<n/2;++i){
    A[i]=A1[i]+wk*A2[i];
    A[i+n/2]=A1[i]-wk*A2[i];
    wk=wk*w1;  
  }
}
int main(){
  scanf("%s%s", s1, s2);
  int n=strlen(s1)-1, m=strlen(s2)-1;
  for(int i=0; i<=n; i++)A[i].x=s1[n-i]-'0';
  for(int i=0; i<=m; i++)B[i].x=s2[m-i]-'0';
  for(m=n+m,n=1;n<=m;n<<=1);
  FFT(A,n,1); FFT(B,n,1);
  for(int i=0;i<n;++i)A[i]=A[i]*B[i];
  FFT(A,n,-1);
  int k=0;
  for(int i=0, t=0; i<n||t; i++){
    t += A[i].x/n+0.5;
    ans[k++] = t%10;
    t /= 10;
  }
  while(k>1 && !ans[k-1]) k--;
  for(int i=k-1; i>=0; i--)printf("%d", ans[i]);
}

 

// 迭代版 4.5s
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
const int N=3e6;
const double PI=acos(-1);

struct complex{
  double x, y;
  complex operator+(const complex& t)const{
    return {x+t.x, y+t.y};}
  complex operator-(const complex& t)const{
    return {x-t.x, y-t.y};}
  complex operator*(const complex& t)const{
    return {x*t.x-y*t.y, x*t.y+y*t.x};}
}A[N], B[N];
char s1[N], s2[N];
int R[N], ans[N];

void FFT(complex A[],int n,int op){
  for(int i=0; i<n; ++i)
    R[i] = R[i/2]/2 + ((i&1)?n/2:0);
  for(int i=0; i<n; ++i)
    if(i<R[i]) swap(A[i],A[R[i]]);
  for(int i=2; i<=n; i<<=1){
    complex w1({cos(2*PI/i),sin(2*PI/i)*op});
    for(int j=0; j<n; j+=i){
      complex wk({1,0});
      for(int k=j; k<j+i/2; ++k){
        complex x=A[k], y=A[k+i/2]*wk;
        A[k]=x+y; A[k+i/2]=x-y; 
        wk=wk*w1;
      }
    }
  }
}
int main(){
  scanf("%s%s", s1, s2);
  int n=strlen(s1)-1, m=strlen(s2)-1;
  for(int i=0; i<=n; i++)A[i].x=s1[n-i]-'0';
  for(int i=0; i<=m; i++)B[i].x=s2[m-i]-'0';
  for(m=n+m,n=1;n<=m;n<<=1);
  FFT(A,n,1); FFT(B,n,1);
  for(int i=0;i<n;++i)A[i]=A[i]*B[i];
  FFT(A,n,-1);
  
  int k=0;
  for(int i=0, t=0; i<n||t; i++){
    t += A[i].x/n+0.5;
    ans[k++] = t%10;
    t /= 10;
  }
  while(k>1 && !ans[k-1]) k--;
  for(int i=k-1; i>=0; i--)printf("%d", ans[i]);
}