#include #include #include #include #include #include #include #include #include #include #include #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 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(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(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 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(xp[pos] >> offset); offset += smallwidth; if(pos + 1 < xn && offset > WordBits) num |= static_cast(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(((long long)xn * stride - 1) / WordBits + 1); Size tmpzn = tmpxn + tmpxn; std::unique_ptr 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(((long long)xn * stride - 1) / WordBits + 1); Size tmpyn = static_cast(((long long)yn * stride - 1) / WordBits + 1); Size tmpzn = tmpxn + tmpyn; std::unique_ptr 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 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 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 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 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 doDP(const int a[6], int n) { int maxX = a[5] * n; vector > dp(n + 1, vector(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 computeLinearRecurrentSequenceCoefficients(const vector &c, long long N) { using namespace mod_poly; Size n = (Size)c.size(); std::unique_ptr 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 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 cntP, cntC; cntP = doDP(primes, P); cntC = doDP(composites, C); int X = P * primes[5] + C * composites[5]; vector ways(X + 1); mod_poly::multiplyPolynomial(ways.data(), X + 1, cntP.data(), (int)cntP.size(), cntC.data(), (int)cntC.size()); vector c(X); for(int i = 0; i < X; ++ i) c[i] = ways[X - i]; vector 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; }