HDU 1402 大数乘法 FFT

大整数乘法，给定两个长度不超过10000的整数，返回乘法的结果。

char*multi(char* number_a, char* number_b)

有疑问欢迎提问，本人学通信的，，手上有 数字信号处理 书，，可以把图搬出来解答的

#include "stdafx.h"
#include <iostream>
#include <cstring>
#include <cmath>
#define pi acos(-1.0)
using namespace std;

struct complex {
	double r, i;
	complex(double real = 0, double image = 0) {
		r = real; i = image;
	}
	complex operator +(const complex &a) {
		return complex(r + a.r, i + a.i);
	}
	complex operator -(const complex &a) {
		return complex(r - a.r, i - a.i);
	}
	complex operator *(const complex &a) {
		return complex(r*a.r - i*a.i, r*a.i + i*a.r);
	}
};

void transAddr(complex *x, int len) {// 变址
	int i, j, k;
	complex temp;
	for (i = 1, j = 0; i < len - 1; i++) {// 0和N-1码位反转结果不变无需交换，略去对其操作。
		// 每次循环把j的开头遇到的1全变为0，随后遇到的第一个0变为1，就能得i的二进制码位反转结果。例如1100-->0010
		k = len >> 1;
		while (j & k) {// 从最高位依次往低位检测，若为1，则通过相与变0
			j &= ~k;
			k >>= 1;
		}
		j |= k;// 第一次遇到的0变为1

		// 交换互为下标码位反转的元素
		if (i < j) {
			temp = x[i];
			x[i] = x[j];
			x[j] = temp;
		}
	}
}

void fft(complex *x, int len, int factor) {// factor == 1时为DFT，factor==-1为IDFT
	int i, j, k;
	complex g, h; // 蝶形运算中的两个零时存储量G(k)、H(k)
	transAddr(x, len); // 先变址
	for (i = 2; i <= len; i <<= 1)// 记多级蝶形运算的每列中各组序列长度为i，从长度为2的序列开始，后一列序列长度为前一列2倍，直到i为原始长度len为止
	{
		complex Wn(cos(-factor * 2 * pi / i), sin(-factor * 2 * pi / i));// 为旋转增量Wn == e^(-2pi/N)按照欧拉公式展开
		for (j = 0; j<len; j += i) {// 某列中各组起始下标
			complex w(1, 0); // 初始化蝴蝶因子。注：Wn^0 == 1
			for (k = j; k < j + i / 2; k++) {// 蝶形运算
				g = x[k];
				h = w * x[k + i / 2];
				x[k] = g + h;
				x[k + i / 2] = g - h;
				w = w * Wn; // 旋转蝴蝶因子
			}
		}
	}
	if (factor == -1) // IDFT跟DFT区别为定义式中蝴蝶因子指数为正 Wn^(-1) == e^(2pi/N)，且大小要乘1/N。注：此处虚部不会再用到，故虚部运算略去。
		for (i = 0; i < len; i++)
			x[i].r /= len;
}

char* multi(char* number_a, char* number_b) {
	// 检查空指针
	if (!number_a || !number_b) {
		cout << "error:输入了空指针\n";
		system("pause");
		return NULL;
	}
	int len_a = strlen(number_a);
	int len_b = strlen(number_b);

	// 检查字符串长度不为空
	if (!len_a || !len_b) {
		cout << "error:输入了空字符串\n";
		system("pause");
		return NULL;
	}

	// 乘法结果的长度最多为2N，且FFT要求len为2的幂
	int len_final = 1;
	while (len_final < 2 * len_a || len_final < 2 * len_b) 
		len_final <<= 1;

	complex *cmplx_a = (complex*)malloc(len_final * sizeof(complex));
	complex *cmplx_b = (complex*)malloc(len_final * sizeof(complex));
	int *product = (int*)malloc(len_final * sizeof(int));
	
	int i,temp;
	for (i = 0; i<len_a; i++) {
		temp = number_a[len_a - 1 - i] - '0';
		if (0 <= temp || temp <= 9) {// 检查输入是否全为数字
			cmplx_a[i] = complex(temp, 0);
		}
		else {
			free(cmplx_a); free(cmplx_b); free(product);
			cout << "error:输入字符串必须全为数字\n";
			system("pause");
			return NULL;
		}
	}
	for (; i < len_final; i++) {
		cmplx_a[i] = complex(0,0);
	}

	for (i = 0; i<len_b; i++) {
		temp = number_b[len_b - 1 - i] - '0';
		if (0 <= temp || temp <= 9) {// 检查输入是否全为数字
			cmplx_b[i] = complex(temp, 0);
		}
		else {
			free(cmplx_a); free(cmplx_b); free(product);
			cout << "error:输入字符串必须全为数字\n";
			system("pause");
			return NULL;
		}
	}
	for (; i < len_final; i++) {
		cmplx_b[i] = complex(0,0);
	}

	fft(cmplx_a, len_final, 1);//通过DFT化卷积为相乘
	fft(cmplx_b, len_final, 1);
	for (i = 0; i < len_final; i++) cmplx_a[i] = cmplx_a[i] * cmplx_b[i];
	fft(cmplx_a, len_final, -1);//逆变换得结果

	// 计算中产生些许误差，故需要四舍五入
	for (i = 0; i < len_final; i++) {
		product[i] = floor(cmplx_a[i].r + 0.5);
	}
  
	// 每个位上的书不超过9，超出则进位
	for (i = 0; i < len_final - 1; i++) {
		product[i + 1] += product[i] / 10;
		product[i] %= 10;
	}

	// 找到最高位在i处
	for (i = len_final - 1; i >= 0; i--) { 
		if (product[i] != 0)
			break;
	}

	// 有i+1位数，再加上\0，共需i+2个单元
	char *number_c = (char*)malloc((i + 2) * sizeof(char));
	char * chrptr = number_c;
	for (; i >= 0; i--) {
		*chrptr++ = product[i] + '0';
	}
	*chrptr = '\0';

	free(cmplx_a); free(cmplx_b); free(product);
	return number_c;
}

int _tmain(int argc, _TCHAR* argv[]) {
	char *number_a = "9999";
	char *number_b = "9999";
	char *number_c = multi(number_a, number_b);
	cout << number_c <<endl;

	system("pause");
	return 0;
}