#include #include #include #include #include #include #include #include #include #include #include "assert.h" using namespace std; #define __MOD__ 1000000007 long long gcd(long long a, long long b){ if(b==0) return a; return gcd(b, a%b); } 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; } /* FFTで使う回転子をある MOD上の原始根にすることで、整数での畳み込みを可能にする O(N log N) まちがい)高速に畳み込むには 原始根^2の冪 = 1 (MOD P) とする必要がある せいかい)高速に畳み込むには、回転子を MOD P 上で位数が 2の冪 となるような x とする必要がある(原始根である必要はない)。 また、この位数はNTTする列より長くなくてはならない。(あらかじめresizeしておく) 実用的には O(10^5) <= 2^k <= O(10^6) となる 17 <= k <= 21 (131072 ~ 2097152) 辺りか 逆変換のときには 回転子^-1 (MOD P) を回転子として使えば良い これを複数の MOD P_k について行い、畳み込んだ後に中国剰余定理を使えば lcm( mod P_k ) 未満の数に復元することができる NTT( f[x] (mod P_k) ) -> F[x] NTT( g[x] (mod P_k) ) -> G[x] Inverse_NTT( F[x]*G[x] (mod P_k) ) -> f*g (mod P_k) Chinese_Remainder_Theorem ( f*g [x] (mod P_k) , MOD Q ) -> f*g [x] (mod Q) */ template 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>= 1){ rbit |= tmp_bit & 1; } rbit >>= 1; if(bit < rbit){ swap(A[bit], A[rbit]); } } int dist = 1; vector& theta = (inverse?inv_theta_v:theta_v); T t = k-1; T half = theta[k-1]; //半回転 for(int level = 1; level= 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= mod) A[i] %= mod; } } return A; } public: const T mod; const T rotater; const T inv_rotater; const T k; vector theta_v; vector 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(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(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 C | C[i] = ΣA[i]B[size-1-j] vector combolution(vector &A, vector &B){ int n = A.size(); assert(A.size() == B.size()); assert( n == (1<= mod) A_NTT[i] %= mod; } return INTT(A_NTT); } }; // 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 x, vector y, long long mod){ int N = x.size(); bool valid = true; //前処理 //gcd(Yi,Yj) == 1 (i!=j) でなくてはならないので、 //共通の因数 g = gcd(Yi,Yj) を見つけたら片側に寄せてしまう /* for(int i=0; i z(N); for(int i=0; i= mod) ans %= mod; tmp *= y[i]; if(tmp >= mod) tmp %= mod; } return ans; } template class kitamasa{ public: int max_n; vector A; vector X; const long long len = 11; //2^18の長さの列について畳み込む vector my_mod = { 98609153,98062337,97792001 }; //NTTする剰余環 vector my_rot = { 1297,29803,179 }; //その単位回転子 vector> NTT; kitamasa(vector& a, vector& x){ for(int i=0; i(my_mod[i],my_rot[i],len)); } //ΣΣ Pi Qj X(i+j) -> Σ Ri Xi ( i = 0...d ) vector func(vector& P, vector& Q){ int d = A.size(); vector tmp(2*d-1, 0); //畳み込む P.resize(1<> R; for(int i=0; i val; for(int j=0; j=d; --i){ for(int j=0; j= __MOD__) tmp[i - d + j] %= __MOD__; } } tmp.resize(d); return tmp; } T calc(T k){ int d = A.size(); int lg = 0; while((1LL<> B(lg, vector(d, 0)); B[0][1] = 1; for(int i=0; i+1 res = B[0]; int r = 0; while(k){ if(k & 1){ res = func(res, B[r]); } k>>=1; r++; } T ret = 0; for(int i=0; i> n >> p >> c; int a[] = {2,3,5,7,11,13}; int b[] = {4,6,8,9,10,12}; int sz = 13*p + 12*c + 1; vector> dp(p+c+1, vector(sz, 0)); dp[0][0] = 1; for(int i=0; i<6; i++){ for(int j=0; j= sz) continue; dp[j+1][k+a[i]] += dp[j][k]; if(dp[j+1][k+a[i]] >= __MOD__) dp[j+1][k+a[i]] %= __MOD__; } } } for(int i=0; i<6; i++){ for(int j=p; j= sz) continue; dp[j+1][k+b[i]] += dp[j][k]; if(dp[j+1][k+b[i]] >= __MOD__) dp[j+1][k+b[i]] %= __MOD__; } } } vector A(sz, 0); for(int i=1; i X(sz, 1); kitamasa ktms(A,X); long long ans = ktms.calc(n+sz-2); cout << ans << endl; return 0; }