FFT/NTT模板既 HDU1402 A * B Problem Plus

@(学习笔记)[FFT, NTT]

Problem Description

Calculate A * B.

Input

Each line will contain two integers A and B. Process to end of file.

Note: the length of each integer will not exceed 50000.

Output

For each case, output A * B in one line.

Sample Input

Sample Output

2

2000

Solution

FFT和NTT都是可以的.

FFT: 没有特别要求.

const int N = 1 << 17;

char str1[N], str2[N];

#include<cstdio>

#include<cstring>

int len1, len2;

struct complex

{

	double real, imag;

	inline complex(){}

	inline complex(double _real, double _imag)

	{

		real = _real, imag = _imag;

	}

	inline friend complex operator *(complex a, complex b)

	{

		return complex(a.real * b.real - a.imag * b.imag, a.imag * b.real + b.imag * a.real);

	}

	inline friend complex operator +(complex a, complex b)

	{

		return complex(a.real + b.real, a.imag + b.imag);

	}

	inline friend complex operator -(complex a, complex b)

	{

		return complex(a.real - b.real, a.imag - b.imag);

	}

}a[N << 1], b[N << 1];

#include<algorithm>

int len;

int rev[N << 1];

inline void prepare(complex *a, complex *b)

{

	int mx = std::max(len1, len2) << 1;

	len = 1;

	int bit = 0;

	while(len < mx)

		len <<= 1, ++ bit;

	rev[0] = 0;

	for(int i = 1; i < len; ++ i)

		rev[i] = rev[i >> 1] >> 1 | (i & 1) << (bit - 1);

}

inline void reverse(complex *a)

{

	for(int i = 0; i < len; ++ i)

		if(rev[i] < i)

			std::swap(a[i], a[rev[i]]);

}

#include<cmath>

const double PI = acos(-1.0);

inline void fft(complex *a, int opt)

{

	reverse(a);

	for(int i = 2; i <= len; i <<= 1)

	{

		complex omega_n = complex(cos(2.0 * PI * (double)opt / (double)i), sin(2.0 * PI * (double)opt / (double)i));

		for(int j = 0; j < len; j += i)

		{

			complex omega = complex(1.0, 0.0); //e ^ 0

			for(int k = j; k < j + i / 2; ++ k)

			{

				complex x = a[k], y = omega * a[k + i / 2];

				a[k] = x + y, a[k + i / 2] = x - y;

				omega = omega * omega_n;

			}

		}

	}

	if(opt == -1)

		for(int i = 0; i < len; ++ i)

			a[i].real /= (double)len;

}

int ans[N << 1];

inline void convolute(complex *a, complex *b)

{

	prepare(a, b);

	fft(a, 1), fft(b, 1);

	for(int i = 0; i < len; ++ i)

		a[i] = a[i] * b[i];

	fft(a, -1);

	for(int i = 0; i < len; ++ i)

		ans[i] = (int)(a[i].real + 0.5);

}

int main()

{

	#ifndef ONLINE_JUDGE

	freopen("HDU1402.in", "r", stdin);

	freopen("HDU1402.out", "w", stdout);

	#endif

	while(gets(str1) && gets(str2))

	{

		len1 = strlen(str1), len2 = strlen(str2);

		memset(a, 0, sizeof(a));

		memset(b, 0, sizeof(b));

		for(int i = 0; i < len1; ++ i)

			a[i] = complex(str1[len1 - i - 1] - '0', 0);

		for(int i = 0; i < len2; ++ i)

			b[i] = complex(str2[len2 - i - 1] - '0', 0);

		convolute(a, b);

		for(int i = 0; i < len; ++ i)

			ans[i + 1] += ans[i] / 10, ans[i] %= 10;

		for(; len && ! ans[len - 1]; -- len);

		for(int i = len - 1; ~ i; -- i)

			putchar(ans[i] + '0');

		if(! len)

			putchar('0');

		putchar('\n');

	}

}

NTT: 要求取模的数\(P\)是费马素数. 对于原根\(g\), 要确保有 \(g ^ 0 \ne g^1 \ne ... \ne g^{p - 2} (modP)\). NTT的过程相当于用 \(\omega_n = g^{ \frac{p - 1}{n}}\) 替代了 \(\omega_n = e^{\frac{2 \cdot \pi}{n} \cdot i}\)的FFT. \(\omega_i\)和\(inv(\omega_i)\)可以在预处理中求出.

const long long P = (479 << 21) + 1;

const long long G = 3;

inline long long quickPower(long long a, long long k, long long mod)

{

	if(! k)

		return 1;

	long long ret = quickPower(a, k >> 1, mod);

	ret = ret * ret % mod;

	if(k & 1)

		ret = ret * a % mod;

	return ret;

}

const long long QUAN = 1 << 5;

long long omega_[QUAN];

inline void getOmega_n()

{

	for(long long i = 0; i < QUAN; ++ i)

		omega_[i] = quickPower(G, (P - 1) / (1 << i), P);

}

const long long LEN = 1 << 16;

char str1[LEN], str2[LEN];

#include<cstdio>

#include<cstring>

long long len1, len2;

long long a[LEN << 1], b[LEN << 1];

#include<algorithm>

long long len;

long long rev[LEN << 1];

inline void prepare()

{

	long long mx = std::max(len1, len2) << 1;

	len = 1;

	long long bit = 0;

	while(len < mx)

		len <<= 1, ++ bit;

	rev[0] = 0;

	for(long long i = 0; i < len; ++ i)

		rev[i] = rev[i >> 1] >> 1 | (i & 1) << (bit - 1);

}

inline void reverse(long long *a)

{

	for(long long i = 0; i < len; ++ i)

		if(rev[i] < i)

			std::swap(a[i], a[rev[i]]);

}

inline void NTT(long long *a, long long opt)

{

	reverse(a);

	long long temp = 0;

	for(long long i = 2; i <= len; i <<= 1)

	{

		++ temp;

		int omega_i = ~ opt ? omega_[temp] : quickPower(omega_[temp], P - 2, P);

		for(long long j = 0; j < len; j += i)

		{

			long long omega = 1;

			for(long long k = j; k < j + i / 2; ++ k)

			{

				long long u = a[k], t = omega * a[k + i / 2] % P;

				a[k] = (u + t) % P, a[k + i / 2] = (u - t + P) % P;

				omega = omega * omega_i % P;

			}

		}

	}

	if(opt == -1)

	{

		long long inv = quickPower(len, P - 2, P);

		for(long long i = 0; i < len; ++ i)

			a[i] = a[i] * inv % P;

	}

}

long long ans[LEN << 1];

inline void convolute(long long *a, long long *b)

{

	prepare();

	NTT(a, 1), NTT(b, 1);

	for(long long i = 0; i < len; ++ i)

		a[i] = a[i] * b[i] % P;

	NTT(a, -1);

	for(int i = 0; i < len; ++ i)

		ans[i] = a[i];

}

int main()

{

	#ifndef ONLINE_JUDGE

	freopen("HDU1402.in", "r", stdin);

	#endif

	getOmega_n();

	while(gets(str1) && gets(str2))

	{

		len1 = strlen(str1), len2 = strlen(str2);

		memset(a, 0, sizeof(a)), memset(b, 0, sizeof(b));

		for(long long i = 0; i < len1; ++ i)

			a[len1 - i - 1] = str1[i] - '0';

		for(long long i = 0; i < len2; ++ i)

			b[len2 - i - 1] = str2[i] - '0';

		convolute(a, b);

		for(int i = 0; i < len; ++ i)

			ans[i + 1] += ans[i] / 10, ans[i] %= 10;

		for(; len && ! ans[len - 1]; -- len);

		for(int i = len - 1; ~ i; -- i)

			putchar(ans[i] + '0');

		if(! len)

			putchar('0');

		putchar('\n');

	}

}

巴特西

FFT/NTT模板既 HDU1402 A * B Problem Plus

Problem Description

Input

Output

Sample Input

Sample Output

Solution

最新文章

热门文章

巴特西

FFT/NTT模板 既 HDU1402 A * B Problem Plus

Problem Description

Input

Output

Sample Input

Sample Output

Solution

最新文章

热门文章

FFT/NTT模板既 HDU1402 A * B Problem Plus