您好,欢迎来到意榕旅游网。
搜索
您的当前位置:首页基于FFT的大整数类高精度乘法

基于FFT的大整数类高精度乘法

来源:意榕旅游网

刘汝佳的《算法竞赛入门经典》中给出过大整数类的代码,但其只实现了输入输出、加减法和比较的功能,并未实现乘法。本文在他的代码的基础上,采用FFT算法实现了高精度乘法运算。
FFT算法讲解:



不熟悉FFT算法的话建议同时浏览对照学习。第一个博客在文章开头所分享的博客给出了较为详细的数学证明,后两个博客对FFT算法做了较为全面的总结,各有千秋。
C++代码如下:

#include <iostream>
#include <algorithm>
#include <complex>
#include <string>
#include <vector>
#define cp complex<double>

const double pi = 3.1415926;
const int maxn = 257;

using namespace std;

class BigInteger
{
private:
	static const int BASE = 100000000;
	static const int WIDTH = 8;

	void recursion_fft(cp *a, int n, int inv) const
	{
		if (n == 1)
			return;
		cp *b = new cp[n];
		int mid = n / 2;
		for (int i = 0; i < mid; i++)
		{
			b[i] = a[i * 2];
			b[i + mid] = a[i * 2 + 1];
		}
		for (int i = 0; i < n; i++)
			a[i] = b[i];
		recursion_fft(a, mid, inv);
		recursion_fft(a + mid, mid, inv);
		for (int i = 0; i < mid; i++)
		{
			cp x(cos(2 * pi * i / n), inv * sin(2 * pi * i / n));
			b[i] = a[i] + x * a[i + mid];
			b[i + mid] = a[i] - x * a[i + mid];
		}
		for (int i = 0; i < n; i++)
			a[i] = b[i];
		delete[] b;
	}

	void iteration_fft(cp *a, int n, int inv) const
	{
		int bit = 0;
		int rev[maxn];
		for (int i = 0; i < n; i++)
			rev[i] = i;
		while ((1 << bit) < n)
			bit++;
		for (int i = 0; i < n; i++)
		{
			rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
			if (i < rev[i])
				swap(a[i], a[rev[i]]);//不加这条if会交换两次(就是没交换)
		}
		for (int mid = 1; mid < n; mid *= 2)//mid是准备合并序列的长度的二分之一
		{
			cp temp(cos(pi / mid), inv * sin(pi / mid));//单位根,pi的系数2已经约掉了
			for (int i = 0; i < n; i += mid * 2)//mid*2是准备合并序列的长度,i是合并到了哪一位
			{
				cp omega(1, 0);
				for (int j = 0; j < mid; j++, omega *= temp)//只扫左半部分,得到右半部分的答案
				{
					cp x = a[i + j], y = omega * a[i + j + mid];
					a[i + j] = x + y, a[i + j + mid] = x - y;//蝴蝶变换
				}
			}
		}
	}

public:
	vector<int> s;

	BigInteger(int num = 0)
	{
		*this = num;
	}

	BigInteger operator = (int num)
	{
		s.clear();
		if (num == 0)
			s.push_back(0);
		while (num > 0)
		{
			s.push_back(num % BASE);
			num /= BASE;
		}
		return *this;
	}

	BigInteger operator = (const string& str)
	{
		s.clear();
		int x, len = (str.length() - 1) / WIDTH + 1;
		for (int i = 0; i < len; i++)
		{
			int last = str.length() - i * WIDTH;
			int start = max(0, last - WIDTH);
			sscanf(str.substr(start, last - start).c_str(), "%d", &x);
			s.push_back(x);
		}
		return *this;
	}

	BigInteger operator + (const BigInteger& b) const
	{
		BigInteger c;
		c.s.clear();
		for (unsigned i = 0, g = 0; ; i++)
		{
			if (g == 0 && i >= s.size() && i >= b.s.size())
				break;
			int x = g;
			if (i < s.size())
				x += s[i];
			if (i < b.s.size())
				x += b.s[i];
			c.s.push_back(x % BASE);
			g = x / BASE;
		}
		return c;
	}

	BigInteger operator * (const BigInteger& bi) const
	{
		string init_a, init_b;
		init_a = BI_to_string();
		init_b = bi.BI_to_string();
		int len_a = init_a.size();
		int len_b = init_b.size();
		int n = 1;
		while (n < len_a + len_b)
			n <<= 1;
		int a[maxn] = { 0 }, b[maxn] = { 0 };
		for (int i = 0; i < len_a; i++)
			a[i] = init_a[len_a - i - 1] - '0';
		for (int i = 0; i < len_b; i++)
			b[i] = init_b[len_b - i - 1] - '0';
		cp cpa[maxn], cpb[maxn];
		for (int i = 0; i < n; i++)
		{
			cpa[i] = cp(a[i], 0);
			cpb[i] = cp(b[i], 0);
		}
		iteration_fft(cpa, n, 1);
		iteration_fft(cpb, n, 1);
		for (int i = 0; i < n; i++)
			cpa[i] *= cpb[i];
		iteration_fft(cpa, n, -1);
		BigInteger c[maxn];
		for (int i = 0; i < n; i++)
			c[i] = (int)(cpa[i].real() / n + 0.5);
		BigInteger ans = 0;
		for (int i = 0; i < n; i++)
			ans = ans + (c[i] << i);
		return ans;
	}

	string BI_to_string() const
	{
		string st;
		st.clear();
		st += to_string(s.back());
		for (int i = s.size() - 2; i >= 0; i--)
		{
			char buf[20];
			sprintf(buf, "%08d", s[i]);
			for (unsigned j = 0; j < strlen(buf); j++)
				st += buf[j];
		}
		return st;
	}

	BigInteger operator << (int num) const //十进制下的左移
	{
		BigInteger x;
		string st = BI_to_string();
		if (st.size() != 1 || st[0] != '0')
			for (int i = 0; i < num; i++)
				st += '0';
		x = st;
		return x;
	}
};

istream& operator >> (istream& in, BigInteger& x)
{
	string s;
	in >> s;
	if (in.fail())
		return in;
	x = s;
	return in;
}

ostream& operator << (ostream& out,const BigInteger& x)
{
	out << x.s.back();
    for(int i = x.s.size()-2; i >= 0; i--)
    {
        char buf[20];
        sprintf(buf, "%08d", x.s[i]);
        for(unsigned j=0; j < strlen(buf); j++)
            out << buf[j];
    }
    return out;
}

int main()
{
	BigInteger a, b;
	cin >> a >> b;
	cout << a * b << endl;
	system("pause");
	return 0;
}

因篇幅问题不能全部显示,请点此查看更多更全内容

Copyright © 2019- yrrf.cn 版权所有 赣ICP备2024042794号-2

违法及侵权请联系:TEL:199 1889 7713 E-MAIL:2724546146@qq.com

本站由北京市万商天勤律师事务所王兴未律师提供法律服务