Submission #420660


Source Code Expand

#include <iostream>
#include <vector>
#include <cstdio>
#include <sstream>
#include <map>
#include <string>
#include <algorithm>
#include <queue>
#include <cmath>
#include <set>
#include <complex>
#include "assert.h"
using namespace std;

long long extgcd(long long a, long long b, long long &x, long long &y){
	long long d=a;
	if(b!=0){
		d = extgcd(b, a%b, y, x);
		y -= (a/b) * x;
	}else{
		x = 1;
		y = 0;
	}
	return d;
}

long long mod_inverse(long long a, long long m){
	long long x,y;
	extgcd(a,m,x,y);
	return (m+x%m)%m;
}

long long mod_pow(long long x, long long y, long long MOD){
	if(x==0) return 0;
	long long ret=1LL;
	while(y>0LL){
		if(y&1LL) ret = (ret * x) % MOD;
		x = (x*x) % MOD;
		y >>= 1LL;
	}
	return ret;
}

template<typename T = long long>
class Number_Theoretic_Transform {
	// return the vector of F[t] or f[x] where
	// F[t] = sum of { f[x] * exp(-j * theta * t * x) } in x = 0 .. N-1 (FFT)
	// f(x) = 1/N * sum of { F[t] * exp(j * theta * t * x) } in t = 0 .. N-1 (inverse FFT)
	// where theta = 2*PI / N
	// N == 2^k
	// use the rotater as (primitive-root of mod) ^ t in NTT, which is used as exp(j*theta)^t in FFT

	//事前に計算した 単位回転子 rotater (MOD mod 上での位数が2^kとなる数) を 引数に与える。
	//逆変換のときには rotater^-1 (MOD mod) を rotaterに与える

	vector< T > do_NTT(vector< T > A, bool inverse){
		int N = A.size();

		//bit reverse
		for(int bit=0; bit<N; bit++){
			int rbit = 0;
			for(int i=0, tmp_bit = bit; i<k-1; i++, rbit <<= 1, tmp_bit >>=  1){
				rbit |= tmp_bit & 1;
			}
			rbit >>= 1;
			if(bit < rbit){
				swap(A[bit], A[rbit]);
			}
		}


		int dist = 1;
		vector<T>& theta = (inverse?inv_theta_v:theta_v);

		T t = k-1;
		T half = theta[k-1];	//半回転

		for(int level = 1; level<k; level++){
			T W_n = theta[t];	//rotater ^ theta (MOD mod)
			T W = 1;							//rotater
			for(int x=0; x < (1<<(level-1)); x++){
				for(int y=x; y+dist < N; y += (dist<<1)){
					T tmp = A[y+dist]*W;
					if(tmp >= mod) tmp %= mod;

					A[y+dist] = A[y] + (tmp*half) % mod;
					if(A[y+dist] >= mod) A[y+dist] %= mod;

					A[y] = A[y] + tmp;
					if(A[y] + tmp >= mod) A[y] %= mod;
				}
				W = W*W_n;
				if(W>=mod) W%=mod;
			}
			dist <<= 1;
			t -= 1;
		}

		if(inverse){
			for(int i=0; i<N; i++){
				A[i] = z * A[i];
				if(A[i] >= mod) A[i] %= mod;
			}
		}
		

		return A;
	}

public:
	const T mod;
	const T rotater;
	const T inv_rotater;
	const T k;
	vector<T> theta_v;
	vector<T> inv_theta_v;
	const T z;

	Number_Theoretic_Transform(T mod_, T rotater_, T k_) : 
		mod(mod_), 
		rotater(rotater_), 
		k(k_), 
		inv_rotater(mod_inverse(rotater_, mod)),
		z(mod_inverse(1<<(k-1), mod)) // 1/Nを掛けるところなので N^-1 MOD modを掛けたいところだけど何故か (N/2)^-1 で上手く行く……
	{

		theta_v = vector<T>(k+1,1);
		theta_v[0] = rotater;
		for(int i=1; i<=k; i++){
			theta_v[i] = theta_v[i-1] * theta_v[i-1];
			if(theta_v[i] >= mod) theta_v[i] %= mod;
		}

		inv_theta_v = vector<T>(k+1,1);
		inv_theta_v[0] = inv_rotater;
		for(int i=1; i<=k; i++){
			inv_theta_v[i] = inv_theta_v[i-1] * inv_theta_v[i-1];
			if(inv_theta_v[i] >= mod) inv_theta_v[i] %= mod;
		}
		
		

	};

	vector< T > NTT(vector< T > A){
		return do_NTT(A, false);
	}

	vector< T > INTT(vector< T > A){
		return do_NTT(A, true);
	}

	// vector<double> C | C[i] = ΣA[i]B[size-1-j]
	vector<T> convolution(vector<T> &A, vector<T> &B){
		int n = A.size();

		assert(A.size() == B.size());
		assert( n == (1<<k) );	//Nは2^kである必要がある(事前にresize)

		auto A_NTT = NTT(A);
		auto B_NTT = NTT(B);

		for(int i=0; i<n; i++){
			A_NTT[i] *= B_NTT[i];
			if(A_NTT[i] >= mod) A_NTT[i] %= mod;
		}
		return INTT(A_NTT);
	}
};

long long gcd(long long a, long long b){
	if(b==0) return a;
	return gcd(b, a%b);
}

long long lcm(long long a, long long b){
	if(a<b) swap(a,b);
	if(b==1) return a;
	return a * (b/gcd(a,b));
}

// Z % Yi = Xi であるようなZを求める。Garnerのアルゴリズム O(N^2)
// 参考 http://techtipshoge.blogspot.jp/2015/02/blog-post_15.html
// http://yukicoder.me/problems/448
long long Chinese_Remainder_Theorem_Garner(vector<long long> x, vector<long long> y, long long MOD){
	int N = x.size();
	//Garner's algorithm
	vector<long long> z(N);
	for(int i=0; i<N; i++){
		z[i] = x[i];
		for(int j=0; j<i; j++){
			z[i] = mod_inverse(y[j], y[i]) % y[i] * (z[i] - z[j]) % y[i];

			z[i] = (z[i]+y[i])%y[i];
		}
	}

	long long ans = 0;
	long long tmp = 1;
	for(int i=0; i<N; i++){
		ans = (ans + z[i] * tmp)%MOD;
		tmp = (tmp * y[i])%MOD;
	}

	return ans;
}



template<typename T = double>
class Fast_Fourier_Transform{
	// return the vector of F[t] or f[x] where
	// F[t] = sum of { f[x] * exp(-j * theta * t * x) } in x = 0 .. N-1 (FFT)
	// f(x) = 1/N * sum of { F[t] * exp(j * theta * t * x) } in t = 0 .. N-1 (inverse FFT)
	// where theta = 2*PI / N
	// N == 2^k
	static vector< complex<T> > do_fft(vector< complex<T> > A, bool inverse){
		const T PI = atan2(1, 0) * 2.0;
		const complex<T> j = complex<T>(0,1); // 0 + j*1.0 (j is the imaginary unit)

		int N = A.size();
		int k = 0;	// N = 2^k
		while(N>>k > 0) k++;
		
		for(int bit=0; bit<N; bit++){
			int rbit = 0;
			for(int i=0, tmp_bit = bit; i<k-1; i++, rbit <<= 1, tmp_bit >>=  1){
				rbit |= tmp_bit & 1;
			}
			rbit >>= 1;
			if(bit < rbit){
				swap(A[bit], A[rbit]);
			}
		}

		int dist = 1;
		T theta = -PI;
		if(inverse) theta = -theta;

		for(int level = 1; level<k; level++){
			complex<T> W_n = exp(j * theta);
			complex<T> W(1,0);
			for(int x=0; x < (1<<level-1); x++){
				for(int y=x; y+dist < N; y += (dist<<1)){
					complex<T> tmp = A[ y+dist ] * W;
					A[ y+dist ] = A[ y ] - tmp;
					A[ y ] += tmp;
				}
				W *= W_n;
			}
			dist <<= 1;
			theta *= 0.5;
		}

		if(inverse){
			T k = 1.0/N;
			for(int i=0; i<N; i++){
				A[i] *= k;
			}
		}
		return A;
	}


	template<typename U=T>
	static void vec_resize(vector<U> &A, const U val){
		int n = A.size();
		int k = 1;
		while(n > k) k<<=1;
		A.resize(k, val);
	}

public:
	Fast_Fourier_Transform(){};

	static vector< complex<T> > FFT(vector< complex<T> > A){
		vec_resize<complex<T>>(A, complex<T>(0,0));
		return do_fft(A, false);
	}

	static vector< complex<T> > IFFT(vector< complex<T> > A){
		//vec_resize<complex<T>>(A, complex<T>(0,0));
		return do_fft(A, true);
	}

	static vector< complex<T> > FFT(vector<T> A){
		vec_resize<T>(A, 0);
		vector<complex<T>> B(A.size());
		for(int i=0; i<B.size(); i++){
			B[i] = complex<T>(A[i], 0);
		}

		return do_fft(B, false);
	}

	static vector< complex<T> > FFT(vector<int> A){
		vec_resize<int>(A, T(0));
		vector<complex<T>> B(A.size());
		for(int i=0; i<B.size(); i++){
			B[i] = complex<T>(A[i], 0);
		}
		return do_fft(B, false);
	}

	static vector<long long> round(vector<complex<T>> &&A){
		vector<long long> ret(A.size());
		for(int i=0; i<A.size(); i++){
			ret[i] = A[i].real() + (A[i].real()<0?-0.5:0.5);
		}
		return ret;
	}

	// vector<double> C | C[i] = ΣA[i]B[j]
	static vector<complex<T>> convolution(vector<T> &A, vector<T> &B){
		//reverse(B.begin(), B.end());

		int n = max(A.size(), B.size());
		A.resize(n, 0);
		B.resize(n, 0);

		auto A_FFT = FFT(A);
		auto B_FFT = FFT(B);
		for(int i=0; i<n; i++){
			A_FFT[i] *= B_FFT[i];
		}
		return IFFT(A_FFT);
	}
};

int main(){
	int n;
	cin >> n;
	/*
	vector<long long> a(1<<18, 0), b(1<<18, 0);
	for(int i=1; i<=n; i++){
		cin >> a[i] >> b[i];
	}
	//reverse(b.begin(), b.end());

	vector<Number_Theoretic_Transform<long long>> my_ntt = {
		Number_Theoretic_Transform<long long>(7340033, 79, 18),
		Number_Theoretic_Transform<long long>(5767169, 127, 18),
		Number_Theoretic_Transform<long long>(8650753, 593, 18)
	};

	vector<vector<long long>> v(3);
	for(int i=0; i<3; i++){
		v[i] = my_ntt[i].convolution(a,b);
	}
	for(int i=1; i<=2*n; i++){
		cout << Chinese_Remainder_Theorem_Garner({v[0][i],v[1][i],v[2][i]}, {7340033,5767169,8650753}, 1LL<<59) << endl;
	}
	
	*/


	vector<double> a(1<<19, 0), b(1<<19, 0);
	for(int i=1; i<=n; i++){
		cin >> a[i] >> b[i];
	}


	auto v = Fast_Fourier_Transform<double>::round(Fast_Fourier_Transform<double>::convolution(a,b));
	for(int i=1; i<=2*n; i++){
		cout << v[i] << endl;
	}


	return 0;
}

Submission Info

Submission Time
Task C - 高速フーリエ変換
User koyumeishi
Language C++11 (GCC 4.9.2)
Score 100
Code Size 8612 Byte
Status AC
Exec Time 1109 ms
Memory 45908 KB

Judge Result

Set Name Sample All
Score / Max Score 0 / 0 100 / 100
Status
AC × 1
AC × 33
Set Name Test Cases
Sample 00_sample_01
All 00_sample_01, 01_00_01, 01_01_19, 01_02_31, 01_03_22, 01_04_31, 01_05_40, 01_06_15, 01_07_39, 01_08_28, 01_09_30, 01_10_23, 01_11_33, 01_12_11, 01_13_28, 01_14_41, 01_15_26, 01_16_49, 01_17_34, 01_18_02, 01_19_33, 01_20_29, 02_00_51254, 02_01_82431, 02_02_17056, 02_03_34866, 02_04_6779, 02_05_65534, 02_06_65535, 02_07_65536, 02_08_65537, 02_09_65538, 02_10_100000
Case Name Status Exec Time Memory
00_sample_01 AC 572 ms 45784 KB
01_00_01 AC 566 ms 45852 KB
01_01_19 AC 565 ms 45848 KB
01_02_31 AC 567 ms 45848 KB
01_03_22 AC 572 ms 45844 KB
01_04_31 AC 585 ms 45848 KB
01_05_40 AC 567 ms 45840 KB
01_06_15 AC 587 ms 45840 KB
01_07_39 AC 610 ms 45840 KB
01_08_28 AC 615 ms 45840 KB
01_09_30 AC 587 ms 45856 KB
01_10_23 AC 560 ms 45840 KB
01_11_33 AC 568 ms 45852 KB
01_12_11 AC 574 ms 45844 KB
01_13_28 AC 573 ms 45856 KB
01_14_41 AC 598 ms 45844 KB
01_15_26 AC 561 ms 45908 KB
01_16_49 AC 560 ms 45852 KB
01_17_34 AC 563 ms 45844 KB
01_18_02 AC 560 ms 45852 KB
01_19_33 AC 562 ms 45844 KB
01_20_29 AC 563 ms 45844 KB
02_00_51254 AC 861 ms 45752 KB
02_01_82431 AC 1104 ms 45744 KB
02_02_17056 AC 649 ms 45852 KB
02_03_34866 AC 790 ms 45852 KB
02_04_6779 AC 596 ms 45844 KB
02_05_65534 AC 946 ms 45860 KB
02_06_65535 AC 915 ms 45852 KB
02_07_65536 AC 920 ms 45840 KB
02_08_65537 AC 921 ms 45848 KB
02_09_65538 AC 912 ms 45848 KB
02_10_100000 AC 1109 ms 45848 KB