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

1
2
1000
2

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('
');
	}
}

NTT: 要求取模的数(P)是费马素数. 对于原根(g), 要确保有 (g ^ 0 e g^1 e ... e 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('
');
	}
}