結果

問題 No.214 素数サイコロと合成数サイコロ (3-Medium)
ユーザー koyumeishikoyumeishi
提出日時 2015-05-29 03:05:58
言語 C++11
(gcc 11.4.0)
結果
AC  
実行時間 1,530 ms / 3,000 ms
コード長 9,282 bytes
コンパイル時間 1,041 ms
コンパイル使用メモリ 93,132 KB
実行使用メモリ 5,932 KB
最終ジャッジ日時 2024-07-06 10:28:26
合計ジャッジ時間 6,404 ms
ジャッジサーバーID
（参考情報） 		judge3 / judge4
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 1,011 ms
5,932 KB
testcase_01 AC 1,434 ms
5,808 KB
testcase_02 AC 1,530 ms
5,928 KB
権限があれば一括ダウンロードができます

ソースコード

diff # 

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

/*
    bool all_zero = true;
    for(int i=0; i<N; i++){
        if(x[i] != 0) all_zero = false;
    }
*/
	//前処理
	//gcd(Yi,Yj) == 1 (i!=j) でなくてはならないので、
	//共通の因数 g = gcd(Yi,Yj) を見つけたら片側に寄せてしまう
	/*
	for(int i=0; i<N; i++){
		for(int j=i+1; j<N; j++){
			if(i == j) continue;
			long long g = gcd(y[i], y[j]);

			if( x[i]%g != x[j]%g ) valid = false;	//解が存在しない

			if(g != 1){
				y[i] /= g; y[j] /= g;
				long long g_ = gcd(y[i], g);
				while(g_ != 1){
					y[i] *= g_;
					g /= g_;
					g_ = gcd(y[i], g);
				}
				y[j] *= g;

				x[i] %= y[i];
				x[j] %= y[j];
			}
		}
	}

	if(!valid){
		cerr << -1 << endl;
		return 0;
	}

	*/
	

	//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 = 13;	//2^18の長さの列について畳み込む

	vector<long long> my_mod = {
		49635329,49586177,49520641/*,49119233*/
	};	//NTTする剰余環
	vector<long long> my_rot = {
		20963,16417,2591/*,5381*/
	};	//その単位回転子
	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__);
		}

		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;
}
0