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-->0010k = len >> 1;while (j & k) {// 从最高位依次往低位检测，若为1，则通过相与变0j &= ~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为IDFTint 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 == 1for (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;}