【Fast Fourier Transform】

文章目录


在这里插入图片描述

在这里插入图片描述

cp omega(int n, int k){
    return cp(cos(2 * PI * k / n), sin(2 * PI * k / n));
}
//omega(n,k)=e^(2*PI*k*i/n)   i指代虚部
// a 是系数数组
void fft(cp  *a, int n, bool inv){
    if(n == 1) return;
    static cp buf[N];
    int m = n / 2;
    for(int i = 0; i < m; i++){ //将每一项按照奇偶分为两组
    // 前半段存放偶数位系数,后半段存放奇数位系数
	    buf[i] = a[2 * i]; 
	    buf[i + m] = a[2 * i + 1];
    }
    for(int i = 0; i < n; i++)
	    a[i] = buf[i];
    fft(a, m, inv); //递归处理两个子问题
    fft(a + m, m, inv);
    for(int i = 0; i < m; i++){ //枚举x,计算A(x)
	    cp x = omega(n, i); 
	    if(inv) x = conj(x); 
	    //conj是一个自带的求共轭复数的函数,精度较高。当复数模为1时,共轭复数等于倒数
	    buf[i] = a[i] + x * a[i + m]; //根据之前推出的结论计算
	    buf[i + m] = a[i] - x * a[i + m];
    }
    for(int i = 0; i < n; i++)
	    a[i] = buf[i];
}

#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <complex>
#define space putchar(' ')
#define enter putchar('\n')
using namespace std;
typedef long long ll;
template <class T>
void read(T &x) {
	char c;
	bool op = 0;
	while (c = getchar(), c < '0' || c > '9')
		if (c == '-') op = 1;
	x = c - '0';
	while (c = getchar(), c >= '0' && c <= '9')
		x = x * 10 + c - '0';
	if (op) x = -x;
}
template <class T>
void write(T x) {
	if (x < 0) putchar('-'), x = -x;
	if (x >= 10) write(x / 10);
	putchar('0' + x % 10);
}
const int N = 1000005;
const double PI = acos(-1);
typedef complex <double> cp; //以后cp就代表 complex<double> 类型
char sa[N], sb[N];
int n = 1, lena, lenb, res[N];
cp a[N], b[N], omg[N], inv[N];
void init() {
	for (int i = 0; i < n; i++) {
		omg[i] = cp(cos(2 * PI * i / n), sin(2 * PI * i / n)); //e^(2*PI*i*i/n) 第一个i指代下标,第二个i指代虚部
		inv[i] = conj(omg[i]);
	}
}
void fft(cp *a, cp *omg) {
	int lim = 0;
	while ((1 << lim) < n) lim++;
	for (int i = 0; i < n; i++) {
		int t = 0;
		for (int j = 0; j < lim; j++)
			if ((i >> j) & 1) t |= (1 << (lim - j - 1));
		if (i < t) swap(a[i], a[t]); // i < t 的限制使得每对点只被交换一次(否则交换两次相当于没交换)
	}
	for (int l = 2; l <= n; l *= 2) {
		int m = l / 2;
		for (cp *p = a; p != a + n; p += l)
			for (int i = 0; i < m; i++) {
				cp t = omg[n / l * i] * p[i + m];
				p[i + m] = p[i] - t;
				p[i] += t;
			}
	}
}
int main() {
	scanf("%s%s", sa, sb);
	lena = strlen(sa), lenb = strlen(sb);
	while (n < lena + lenb) n *= 2;
	for (int i = 0; i < lena; i++)
		a[i].real(sa[lena - 1 - i] - '0');  //将字符串倒序放入a数组中
	for (int i = 0; i < lenb; i++)
		b[i].real(sb[lenb - 1 - i] - '0');   //将字符串倒序放入b数组中
	init();
	fft(a, omg);
	fft(b, omg);
	for (int i = 0; i < n; i++)
		a[i] *= b[i];
	fft(a, inv);
	for (int i = 0; i < n; i++) {
		res[i] += floor(a[i].real() / n + 0.5);
		res[i + 1] += res[i] / 10;
		res[i] %= 10;
	}
	for (int i = res[lena + lenb - 1] ? lena + lenb - 1 : lena + lenb - 2; i >= 0; i--)
		putchar('0' + res[i]);
	enter;
	system("pause");
	return 0;
}

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值