結果
| 問題 | No.215 素数サイコロと合成数サイコロ (3-Hard) |
| ユーザー |
 koyumeishi
|
| 提出日時 | 2015-05-29 03:54:50 |
| 言語 | C++11 (gcc 11.4.0) |
| 結果 |
TLE
|
| 実行時間 | - |
| コード長 | 8,863 bytes |
| コンパイル時間 | 1,851 ms |
| コンパイル使用メモリ | 94,608 KB |
| 実行使用メモリ | 48,304 KB |
| 最終ジャッジ日時 | 2023-09-20 15:26:41 |
| 合計ジャッジ時間 | 14,349 ms |
|
ジャッジサーバーID （参考情報） |
judge14 / judge11 |
（要ログイン）
ソースコード
#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;
#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<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> combolution(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);
}
};
// 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();
bool valid = true;
//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 += z[i]*tmp;
if(ans >= mod) ans %= mod;
tmp *= y[i];
if(tmp >= mod) tmp %= mod;
}
/*
if(all_zero){
long long lcm_ = 1;
for(int i=0; i<N; i++){
lcm_ = (lcm_ * y[i]%mod)%mod;
}
ans = lcm_;
}
*/
return ans;
}
template<class T>
class kitamasa{
public:
int max_n;
vector<T> A;
vector<T> X;
const long long len = 15; //2^18の長さの列について畳み込む
vector<long long> my_mod = {
98861057,97615873,97320961/*,96010241,95027201,94863361,94273537,93880321,93192193*/
}; //NTTする剰余環
vector<long long> my_rot = {
9859,6427,11393/*,157,6701,6637,17299,2399,8329*/
}; //その単位回転子
vector<Number_Theoretic_Transform<long long>> NTT;
kitamasa(vector<T>& a, vector<T>& x){
for(int i=0; i<a.size(); i++) A.push_back(a[i]);
for(int i=0; i<x.size(); i++) X.push_back(x[i]);
for(int i=0; i<my_mod.size(); i++) NTT.push_back(Number_Theoretic_Transform<long long>(my_mod[i],my_rot[i],len));
}
//ΣΣ Pi Qj X(i+j) -> Σ Ri Xi ( i = 0...d )
vector<T> func(vector<T> P, vector<T> Q){
int d = A.size();
vector<T> tmp(2*d-1, 0);
//畳み込む
P.resize(1<<len);
Q.resize(1<<len);
vector<vector<long long>> R;
for(int i=0; i<my_mod.size(); i++){
R.push_back(NTT[i].combolution(P,Q));
}
for(int i=0; i<tmp.size(); i++){
vector<long long> val;
for(int j=0; j<R.size(); j++){
val.push_back(R[j][i]);
}
tmp[i] = Chinese_Remainder_Theorem_Garner(val, my_mod, __MOD__);
}
//
// ここを O(n) もしくは O(n log n) 程度に落とす必要がある
// エスケープしててもTLEなようだと他も改善する必要あり
//
//ここの O(N^2)をどげんかせんといかん
//for(int i=2*d-2; i>=d; --i){
// if(tmp[i] == 0) continue;
// for(int j=0; j<d; ++j){
// tmp[i - d + j] += (tmp[i] * A[j]);
// if(tmp[i - d + 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<<lg) <= k) lg++;
// X[2^n] = Σ B[n][i] X[i] (i = 0...d)
vector<vector<T>> B(lg, vector<T>(d, 0));
B[0][1] = 1;
for(int i=0; i+1<lg; i++){
B[i+1] = func(B[i], B[i]);
}
vector<T> 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<d; i++){
ret += (res[i] * X[i]);
ret %= __MOD__;
}
return ret;
}
};
int main(){
long long n;
int p,c;
cin >> 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<vector<long long>> dp(p+c+1, vector<long long>(sz, 0));
dp[0][0] = 1;
for(int i=0; i<6; i++){
for(int j=0; j<p; j++){
for(int k=0; k<sz; k++){
if(k+a[i] >= 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<p+c; j++){
for(int k=0; k<sz; k++){
if(k+b[i] >= 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<long long> A(sz, 0);
for(int i=1; i<sz; i++){
A[i-1] = dp[p+c][i];
}
reverse(A.begin(), A.end());
vector<long long> X(sz, 1);
kitamasa<long long> ktms(A,X);
long long ans = ktms.calc(n+sz-2);
cout << ans << endl;
return 0;
}