快速傅里叶变换实现两个超大整数(100000位)乘法

[eiclpy]

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
const int SIZEN = 800010, BASE = 10;
const double PI = acos(-1), eps = 1e-6;
class Complex
{ //复数
  public:
	double real, imag;
	Complex(double x = 0, double y = 0) { real = x, imag = y; }
};
Complex operator+(Complex a, Complex b) { return Complex(a.real + b.real, a.imag + b.imag); }
Complex operator-(Complex a, Complex b) { return Complex(a.real - b.real, a.imag - b.imag); }
Complex operator*(Complex a, Complex b) { return Complex(a.real * b.real - a.imag * b.imag, a.real * b.imag + b.real * a.imag); }
Complex operator/(Complex a, double b) { return Complex(a.real / b, a.imag / b); }
class Poly
{
  public:
	int n; //n-1次多项式,有n个系数
	Complex s[SIZEN];
	void Initialize(char str[])
	{
		n = strlen(str);
		for (int i = 0; i < n; i++)
			s[i] = Complex(str[n - 1 - i] - '0', 0);
	}
	void Assign(char str[])
	{
		static int a[SIZEN];
		int len;
		for (int i = 0; i < n; i++)
			a[i] = int(s[i].real + 0.5);
		for (len = 0; len < n || a[len]; len++)
		{
			a[len + 1] += a[len] / BASE;
			a[len] %= BASE;
		}
		while (len > 1 && !a[len - 1])
			len--;
		for (int i = 0; i < len; i++)
			str[i] = a[len - 1 - i] + '0';
		str[len] = 0;
	}
	void Print(void)
	{
		static char str[SIZEN];
		Assign(str);
		printf("%s\n", str);
	}
	void Rader_Transform(void)
	{
		int j = 0, k;
		for (int i = 0; i < n; i++)
		{
			if (j > i)
				swap(s[i], s[j]);
			k = n;
			//进位
			while (j & (k >>= 1))
				j &= ~k;
			j |= k;
		}
	}
	void FFT(bool type)
	{								 //type=1是DFT求值(系数求点),type=0是IDFT插值(点求系数)
		Rader_Transform();			 //重新排列数组
		double pi = type ? PI : -PI; //IDFT的要点1:π变负
		Complex w0;
		for (int step = 1; step < n; step <<= 1)
		{
			//在这层循环中我们把相邻两个长度为step的点值表达式合并
			for (int H = 0; H < n; H += step << 1)
			{
				//将[H,H+step)和[H+step,H+2*step)这两个点值表达式合为一个
				double alpha = pi / step; //2*step次单位根,转动的角度为alpha的倍数
				Complex wn(cos(alpha), sin(alpha)), wk(1.0, 0.0);
				//wn是幅角为alpha的单位根,wk从1开始每次乘以wn就能遍历每个2*step次单位根
				for (int k = 0; k < step; k++)
				{
					//执行一次旋转操作
					int Ek = H + k, Ok = H + k + step; //旋转操作的第一项和第二项
					Complex t = wk * s[Ok];			   //旋转因子
					s[Ok] = s[Ek] - t;
					s[Ek] = s[Ek] + t;
					wk = wk * wn;
				}
			}
		}
		if (!type)
			for (int i = 0; i < n; i++)
				s[i] = s[i] / n; //IDFT的要点2:除以n
	}
	void operator*=(Poly &b)
	{
		int S = 1;
		while (S < n + b.n)
			S <<= 1;
		n = b.n = S;
		FFT(true);
		b.FFT(true);
		for (int i = 0; i < n; i++)
			s[i] = s[i] * b.s[i];
		FFT(false);
	}
};
Poly A, B;
int main()
{
	char str[SIZEN];
	scanf("%s", str);
	A.Initialize(str);
	scanf("%s", str);
	B.Initialize(str);
	A *= B;
	A.Print();
	return 0;
}

[eiclpy]

测试

import random
import os

for _ in range(10):
    with open('input.txt', 'w') as f:
        val1 = random.randint(9*10**10000, 10**100000)
        val2 = random.randint(9*10**10000, 10**100000)
        print(val1, file=f)
        print(val2, file=f)
    os.system(r'.\fft.exe < input.txt > out.txt')
    with open('out.txt', 'r') as f1:
        if f1.read().strip() == str(val1*val2):
            print('OK')
        else:
            print("Wrong")