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 |
|
|
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 |