#67954 (C++11（廃止可能性あり）) No.215 素数サイコロと合成数サイコロ (3-Hard)

• 問題: No.215 素数サイコロと合成数サイコロ (3-Hard)
• ユーザー: anta
• 提出日時: 2015-12-22 20:08:19
• 言語: C++11（廃止可能性あり） (gcc 13.3.0)
• 結果: CE (最新)
  AC
  (最初)
• 実行時間: -
• コード長: 11,325 bytes
• コンパイル時間: 1,914 ms
• コンパイル使用メモリ: 86,148 KB
• 最終ジャッジ日時: 2024-11-14 19:31:27
• 合計ジャッジ時間: 2,271 ms
• ジャッジサーバーID （参考情報）: judge4 / judge2

コンパイルメッセージ

/usr/bin/ld: /tmp/ccyZ6CGt.o: in function `DynamicLibrary::get(char const*) const [clone .isra.0]':
main.cpp:(.text.startup+0x5): undefined reference to `__libc_dlsym'
/usr/bin/ld: /tmp/ccyZ6CGt.o: in function `_GLOBAL__sub_I_gmp':
main.cpp:(.text.startup+0x457): undefined reference to `__libc_dlopen_mode'
collect2: error: ld returned 1 exit status

ソースコード

#include <string>
#include <iostream>
#include <cstdlib>
#include <tuple>
#include <algorithm>
#include <cstdint>
#include <cassert>
#include <memory>
#include <vector>
#include <numeric>

#include <gmp.h>

#if defined(_WIN32)
extern "C" {
	void* __stdcall LoadLibraryA(const char *name);
	void* __stdcall GetProcAddress(void *handle, const char *name);
}
struct DynamicLibrary {
	DynamicLibrary(const char *name) { _handle = LoadLibraryA(name); }
	void *get(const char *name) const { void *p = GetProcAddress(_handle, name); if(!p) std::cerr << "can't find: " << name << std::endl, std::abort(); return p; }
	void *_handle;
} gmp("mpir.dll");
#else
extern "C" {
	void* __libc_dlopen_mode(const char *x, int y);
	void* __libc_dlsym(void *x, const char *y);
}
struct DynamicLibrary {
	DynamicLibrary(const char *name) { _handle = __libc_dlopen_mode(name, 2); }
	void *get(const char *name) const { void *p = __libc_dlsym(_handle, name); if(!p) std::cerr << "can't find: " << name << std::endl, std::abort(); return p; }
	void *_handle;
} gmp("/usr/lib64/libgmp.so.10");	//yukicoder
#endif

#define EXPAND_MACRO_TO_STR_2(a) #a
#define EXPAND_MACRO_TO_STR(a) EXPAND_MACRO_TO_STR_2(a)
#define DECL_GMP(name) const auto my_##name = (decltype(::name)*)gmp.get(EXPAND_MACRO_TO_STR(name))
DECL_GMP(mpn_mul);
DECL_GMP(mpn_sqr);

template<int MOD>
struct ModInt {
	static const int Mod = MOD;
	unsigned x;
	ModInt() : x(0) {}
	ModInt(signed sig) { int sigt = sig % MOD; if(sigt < 0) sigt += MOD; x = sigt; }
	ModInt(unsigned unsig) { x = unsig % MOD; }
	ModInt(signed long long sig) { int sigt = sig % MOD; if(sigt < 0) sigt += MOD; x = sigt; }
	int get() const { return (int)x; }

	ModInt &operator+=(ModInt that) { if((x += that.x) >= MOD) x -= MOD; return *this; }
	ModInt &operator-=(ModInt that) { if((x += MOD - that.x) >= MOD) x -= MOD; return *this; }
	ModInt &operator*=(ModInt that) { x = (unsigned long long)x * that.x % MOD; return *this; }

	ModInt operator+(ModInt that) const { return ModInt(*this) += that; }
	ModInt operator-(ModInt that) const { return ModInt(*this) -= that; }
	ModInt operator*(ModInt that) const { return ModInt(*this) *= that; }

	bool operator==(ModInt that) const { return x == that.x; }
	bool operator!=(ModInt that) const { return x != that.x; }
	ModInt operator-() const { ModInt t; t.x = x == 0 ? 0 : Mod - x; return t; }

};
typedef ModInt<1000000007> mint;

namespace mod_poly {

inline int getBitSize(unsigned x) {
#ifndef __GNUC__
	unsigned long res;
	_BitScanReverse(&res, x);
	return static_cast<int>(res) + 1;
#else
	return 32 - __builtin_clz(x);
#endif
}

enum { WordBits = sizeof(mp_limb_t) * 8 };
typedef unsigned Num;
enum { NumBits = sizeof(Num) * 8 };

typedef int Size;
typedef const mint *SrcPoly;
typedef mint *DstPoly, SrcDstPoly;
typedef const mp_limb_t *SrcPtr;
typedef mp_limb_t *DstPtr, *SrcDstPtr;

void zeroPolynomial(DstPoly resp, Size resn) {
	for(Size i = 0; i < resn; ++ i)
		resp[i] = mint();
}

void addPolynomial(DstPoly resp, SrcPoly xp, SrcPoly yp, Size xyn) {
	for(Size i = 0; i < xyn; ++ i)
		resp[i] = xp[i] + yp[i];
}
void subtractPolynomial(DstPoly resp, SrcPoly xp, SrcPoly yp, Size xyn) {
	for(Size i = 0; i < xyn; ++ i)
		resp[i] = xp[i] - yp[i];
}

void multiplyConstantPolynomial(DstPoly resp, SrcPoly xp, Size xn, mint y) {
	for(Size i = 0; i < xn; ++ i)
		resp[i] = xp[i] * y;
}

void packToBits(DstPtr resp, int resn, SrcPoly xp, Size xn, int width, int stride) {
	if(width > WordBits) {
		std::cerr << "unimplemented" << std::endl;
		std::abort();
	}
	for(Size i = 0; i < resn; ++ i)
		resp[i] = 0;
	Size pos = 0;
	int sa = stride / WordBits, sb = stride % WordBits;
	int offset = 0;
	for(Size i = 0; i < xn; ++ i) {
		assert(pos < resn);
		mp_limb_t x = static_cast<mp_limb_t>(xp[i].get());
		int right = offset + width;
		resp[pos] |= x << offset;
		if(right > WordBits) {
			assert(pos + 1 < resn);
			resp[pos + 1] |= x >> (WordBits - offset);
		}
		pos += sa;
		if((offset += sb) >= WordBits)
			offset -= WordBits, ++ pos;
	}
}

void unpackFromBitsMod(DstPoly resp, Size resn, SrcPtr xp, Size xn, int stride) {
	assert(stride > 0);
	int K = (stride - 1) / NumBits + 1, rems = stride - (K - 1) * NumBits;
	Num remmask = ((Num)2 << (rems - 1)) - 1;
	const mint Base = mint(~0U) + 1;
	std::unique_ptr<Num[]> buf(new Num[K]);
	Num *nums = buf.get();
	Size pos = 0, offset = 0;
	for(Size i = 0; i < resn; ++ i) {
		for(int k = 0; k < K; ++ k) {
			int smallwidth = k < K - 1 ? NumBits : rems;
			if(pos >= xn) {
				nums[k] = 0;
				offset += smallwidth;
				continue;
			}
			Num num = static_cast<Num>(xp[pos] >> offset);
			offset += smallwidth;
			if(pos + 1 < xn && offset > WordBits)
				num |= static_cast<Num>(xp[pos + 1] << (WordBits - (offset - smallwidth)));
			nums[k] = num;
			if(offset >= WordBits)
				offset -= WordBits, ++ pos;
		}
		nums[K - 1] &= remmask;
		mint cur = mint(nums[K - 1]);
		for(int k = K - 2; k >= 0; -- k) {
			cur *= Base;
			cur += nums[k];
		}
		resp[i] = cur;
	}
}

void multiplyPolynomial(DstPoly resp, Size resn, SrcPoly xp, Size xn, SrcPoly yp, Size yn) {
	if(xn > resn) xn = resn;
	if(yn > resn) yn = resn;
	if(resn > xn + yn - 1) {
		zeroPolynomial(resp + (xn + yn - 1), resn - (xn + yn - 1));
		resn = xn + yn - 1;
	}
	if(xn < yn) std::swap(xp, yp), std::swap(xn, yn);
	if(yn == 0) {
		zeroPolynomial(resp, resn);
		return;
	}
	if(yn == 1) {
		multiplyConstantPolynomial(resp, xp, xn, yp[0]);
		zeroPolynomial(resp + xn, resn - xn);
		return;
	}
	int width = getBitSize(mint::Mod - 1);
	int stride = width * 2 + getBitSize(yn - 1);
	if(xp == yp && xn == yn) {
		Size tmpxn = static_cast<Size>(((long long)xn * stride - 1) / WordBits + 1);
		Size tmpzn = tmpxn + tmpxn;
		std::unique_ptr<mp_limb_t[]> tmpbuf(new mp_limb_t[tmpxn + tmpzn]);
		mp_limb_t *tmpx = tmpbuf.get(), *tmpz = tmpx + tmpxn;
		packToBits(tmpx, tmpxn, xp, xn, width, stride);
		my_mpn_sqr(tmpz, tmpx, tmpxn);
		unpackFromBitsMod(resp, resn, tmpz, tmpzn, stride);
	} else {
		Size tmpxn = static_cast<Size>(((long long)xn * stride - 1) / WordBits + 1);
		Size tmpyn = static_cast<Size>(((long long)yn * stride - 1) / WordBits + 1);
		Size tmpzn = tmpxn + tmpyn;
		std::unique_ptr<mp_limb_t[]> tmpbuf(new mp_limb_t[tmpxn + tmpyn + tmpzn]);
		mp_limb_t *tmpx = tmpbuf.get(), *tmpy = tmpx + tmpxn, *tmpz = tmpy + tmpyn;
		packToBits(tmpx, tmpxn, xp, xn, width, stride);
		packToBits(tmpy, tmpyn, yp, yn, width, stride);
		my_mpn_mul(tmpz, tmpx, tmpxn, tmpy, tmpyn);
		unpackFromBitsMod(resp, resn, tmpz, tmpzn, stride);
	}
}

void squarePolynomial(DstPoly resp, Size resn, SrcPoly xp, Size xn) {
	return multiplyPolynomial(resp, resn, xp, xn, xp, xn);
}


void invertPowerSeries(DstPoly resp, Size resn, SrcPoly xp, Size xn) {
	if(resn == 0) return;
	assert(resp != xp);
	mint two = 2;
	std::unique_ptr<mint[]> tmpbuf(new mint[resn + resn]);
	mint *tmp1 = tmpbuf.get(), *tmp2 = tmp1 + resn;
	zeroPolynomial(resp, resn);
	resp[0] = xp[0];
	Size curn = 1;
	while(curn < resn) {
		Size nextn = std::min(resn, curn * 2);
		squarePolynomial(tmp1, nextn, resp, curn);
		multiplyPolynomial(tmp2, nextn, tmp1, nextn, xp, std::min(nextn, xn));
		multiplyConstantPolynomial(resp, resp, curn, two);
		subtractPolynomial(resp, resp, tmp2, nextn);
		curn = nextn;
	}
}

void reversePolynomial(DstPoly resp, SrcPoly xp, Size xn) {
	if(resp == xp) {
		std::reverse(resp, resp + xn);
	} else {
		for(Size i = 0; i < xn; ++ i)
			resp[xn - 1 - i] = xp[i];
	}
}

void precomputePolynomialInverse(DstPoly resp, Size resn, SrcPoly xp, Size xn) {
	std::unique_ptr<mint[]> tmp(new mint[xn]);
	reversePolynomial(tmp.get(), xp, xn);
	invertPowerSeries(resp, resn, tmp.get(), xn);
}

void dividePolynomialPrecomputed(DstPoly resp, Size resn, SrcPoly revxp, Size xn, SrcPoly invp) {
	multiplyPolynomial(resp, resn, revxp, xn, invp, resn);
	reversePolynomial(resp, resp, resn);
}

void remainderPolynomialPrecomputed(DstPoly resp, SrcPoly xp, Size xn, SrcPoly yp, Size yn, SrcPoly invp) {
	if(xn < yn) {
		for(Size i = 0; i < xn; ++ i)
			resp[i] = xp[i];
		zeroPolynomial(resp + xn, yn - 1 - xn);
		return;
	}
	if(yn == 0)
		abort();
	if(yp[yn - 1].get() != 1) {
		std::cerr << "unimplemented" << std::endl;
		abort();
	}
	if(yn == 1) return;
	Size quotn = xn - yn + 1;
	std::unique_ptr<mint[]> tmpbuf(new mint[xn + quotn + yn - 1]);
	mint *revxp = tmpbuf.get(), *quot = revxp + xn, *quotmul = quot + quotn;
	reversePolynomial(revxp, xp, xn);
	dividePolynomialPrecomputed(quot, quotn, revxp, xn, invp);
	multiplyPolynomial(quotmul, yn - 1, quot, std::min(quotn, yn - 1), yp, yn - 1);
	subtractPolynomial(resp, xp, quotmul, yn - 1);
}

void polynomialPowerModX(DstPoly resp, long long K, SrcPoly yp, Size yn, SrcPoly invp) {
	assert(K >= 0 && yn > 0);
	assert(yp[yn - 1] == 1);
	if(yn == 1) return;
	if(K == 0) {
		resp[0] = 1;
		zeroPolynomial(resp + 1, yn - 2);
		return;
	}
	int l = 0;
	while((K >> l) > 1) ++ l;
	zeroPolynomial(resp, yn - 1);
	resp[1] = 1;
	std::unique_ptr<mint[]> tmpbuf(new mint[yn * 2 - 1]);
	mint *square = tmpbuf.get();
	bool init = true;
	Size curdeg = 1, degy = yn - 1;
	for(-- l; l >= 0; -- l) {
		squarePolynomial(square, curdeg * 2 + 1, resp, curdeg + 1);
		curdeg *= 2;
		remainderPolynomialPrecomputed(resp, square, curdeg + 1, yp, yn, invp);
		if(curdeg >= degy) curdeg = degy - 1;
		if(K >> l & 1) {
			if(curdeg < degy - 1) {
				for(Size i = curdeg + 1; i > 0; -- i)
					resp[i] = resp[i - 1];
				++ curdeg;
			} else {
				mint neglc = -resp[degy - 1];
				for(Size i = degy - 1; i > 0; -- i)
					resp[i] = resp[i - 1] + yp[i] * neglc;
				resp[0] = yp[0] * neglc;
			}
		}
	}
}


}	//namespace mod_poly

using std::vector;
using std::cin;
vector<mint> doDP(const int a[6], int n) {
	int maxX = a[5] * n;
	vector<vector<mint> > dp(n + 1, vector<mint>(maxX + 1));
	dp[0][0] = mint(1);
	int z = a[0];
	for(int k = 0; k < 6; ++ k) {
		int d = a[k];
		for(int t = 0; t < n; ++ t) {
			for(int i = t * z; i <= t * d; ++ i) {
				mint x = dp[t][i];
				if(x.x != 0)
					dp[t + 1][i + d] += x;
			}
		}
	}
	return dp[n];
}

//a_{i+n} = \sum_{j=0}^{n-1} c[j] a_{i+j}
vector<mint> computeLinearRecurrentSequenceCoefficients(const vector<mint> &c, long long N) {
	using namespace mod_poly;
	Size n = (Size)c.size();
	std::unique_ptr<mint[]> tmpbuf(new mint[(n + 1) * 3]);
	mint *poly = tmpbuf.get(), *inv = poly + (n + 1), *resp = inv + (n + 1);
	for(Size i = 0; i < n; ++ i)
		poly[i] = -c[i];
	poly[n] = 1;
	precomputePolynomialInverse(inv, n + 1, poly, n + 1);
	polynomialPowerModX(resp, N, poly, n + 1, inv);
	vector<mint> res(n);
	for(Size i = 0; i < n; ++ i)
		res[i] = resp[i];
	return res;
}

int main() {
	long long N; int P; int C;
	while(cin >> N >> P >> C) {
		const int primes[6] = {2, 3, 5, 7, 11, 13};
		const int composites[6] = {4, 6, 8, 9, 10, 12};
		vector<mint> cntP, cntC;
		cntP = doDP(primes, P);
		cntC = doDP(composites, C);
		int X = P * primes[5] + C * composites[5];
		vector<mint> ways(X + 1);
		mod_poly::multiplyPolynomial(ways.data(), X + 1, cntP.data(), (int)cntP.size(), cntC.data(), (int)cntC.size());
		vector<mint> c(X);
		for(int i = 0; i < X; ++ i)
			c[i] = ways[X - i];
		vector<mint> coeffs;
		coeffs = computeLinearRecurrentSequenceCoefficients(c, X + (N - 1));
		mint ans = std::accumulate(coeffs.begin(), coeffs.end(), mint());
		std::printf("%d\n", (int)ans.get());
	}
	return 0;
}