Submission #420592
Source Code Expand
#include <iostream>
#include <vector>
#include <cstdio>
#include <sstream>
#include <map>
#include <string>
#include <algorithm>
#include <queue>
#include <cmath>
#include <set>
#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;
}
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;
}
return 0;
}
Submission Info
Submission Time |
|
Task |
C - 高速フーリエ変換 |
User |
koyumeishi |
Language |
C++11 (GCC 4.9.2) |
Score |
0 |
Code Size |
5484 Byte |
Status |
WA |
Exec Time |
2398 ms |
Memory |
17276 KB |
Judge Result
Set Name |
Sample |
All |
Score / Max Score |
0 / 0 |
0 / 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 |
952 ms |
17168 KB |
01_00_01 |
AC |
873 ms |
17180 KB |
01_01_19 |
AC |
1017 ms |
17192 KB |
01_02_31 |
AC |
1045 ms |
17180 KB |
01_03_22 |
AC |
1027 ms |
17184 KB |
01_04_31 |
AC |
1047 ms |
17148 KB |
01_05_40 |
AC |
1059 ms |
17180 KB |
01_06_15 |
AC |
1004 ms |
17188 KB |
01_07_39 |
AC |
1059 ms |
17184 KB |
01_08_28 |
AC |
1044 ms |
17276 KB |
01_09_30 |
AC |
1044 ms |
17112 KB |
01_10_23 |
AC |
1031 ms |
17176 KB |
01_11_33 |
AC |
1051 ms |
17176 KB |
01_12_11 |
AC |
993 ms |
17168 KB |
01_13_28 |
AC |
1040 ms |
17176 KB |
01_14_41 |
AC |
1062 ms |
17184 KB |
01_15_26 |
AC |
1034 ms |
17180 KB |
01_16_49 |
AC |
1082 ms |
17144 KB |
01_17_34 |
AC |
1068 ms |
17176 KB |
01_18_02 |
AC |
906 ms |
17188 KB |
01_19_33 |
AC |
1052 ms |
17164 KB |
01_20_29 |
AC |
1042 ms |
17176 KB |
02_00_51254 |
AC |
1949 ms |
17180 KB |
02_01_82431 |
WA |
2170 ms |
17180 KB |
02_02_17056 |
AC |
1535 ms |
17176 KB |
02_03_34866 |
AC |
1739 ms |
17172 KB |
02_04_6779 |
AC |
1415 ms |
17176 KB |
02_05_65534 |
AC |
2103 ms |
17180 KB |
02_06_65535 |
AC |
2087 ms |
17176 KB |
02_07_65536 |
WA |
2094 ms |
17176 KB |
02_08_65537 |
WA |
2103 ms |
17184 KB |
02_09_65538 |
WA |
2090 ms |
17180 KB |
02_10_100000 |
WA |
2398 ms |
17204 KB |