#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifndef __GNUC__ #include #endif #define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i)) #define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i)) #define reu(i,l,u) for(int (i)=(int)(l);(i)<(int)(u);++(i)) #if defined(_MSC_VER) || __cplusplus > 199711L #define aut(r,v) auto r = (v) #else #define aut(r,v) __typeof(v) r = (v) #endif #define each(it,o) for(aut(it, (o).begin()); it != (o).end(); ++ it) #define all(o) (o).begin(), (o).end() #define pb(x) push_back(x) #define mp(x,y) make_pair((x),(y)) #define mset(m,v) memset(m,v,sizeof(m)) #define INF 0x3f3f3f3f #define INFL 0x3f3f3f3f3f3f3f3fLL using namespace std; typedef vector vi; typedef pair pii; typedef vector > vpii; typedef long long ll; template inline void amin(T &x, U y) { if(y < x) x = y; } template inline void amax(T &x, U y) { if(x < y) x = y; } #ifdef MY_LOCAL_RUN #include "C:\Dropbox\backup\implements\Util\MyAssert.hpp" #undef assert #define assert my_assert #endif namespace multiprecision { #if defined(_M_X64) || defined(__amd64__) #define MP_64 #endif #ifdef MP_64 typedef unsigned long long Word; #else typedef unsigned Word; #endif typedef int Size; //max size * WordBits typedef int BitSize; typedef const Word *SrcPtr; typedef Word *DstPtr, *SrcDstPtr; const int WordBits = sizeof(Word) * 8; const Word WordMax = ~static_cast(0); #define srep(i,n) for(Size (i)=0;(i)<(Size)(n);++(i)) inline int getBitSize(Word x) { #ifdef MP_64 #ifndef __GNUC__ unsigned long res; _BitScanReverse64(&res, x); return static_cast(res) + 1; #else return 64 - __builtin_clzll(x); #endif #else #ifndef __GNUC__ unsigned long res; _BitScanReverse(&res, x); return static_cast(res)+ 1; #else return 32 - __builtin_clz(x); #endif #endif } inline void multiplyWords(Word &resLo, Word &resHi, Word x, Word y) { #ifdef MP_64 #ifndef __GNUC__ resLo = _umul128(x, y, &resHi); #else unsigned __int128 z = (unsigned __int128)x * y; resLo = (Word)z; resHi = (Word)(z >> 64); #endif #else unsigned long long z = (unsigned long long)x * y; resLo = (unsigned)z; resHi = (unsigned)(z >> 32); #endif } inline Word multiplyWordsLo(Word x, Word y) { Word lo, hi; multiplyWords(lo, hi, x, y); return lo; } inline Word multiplyWordsHi(Word x, Word y) { Word lo, hi; multiplyWords(lo, hi, x, y); return hi; } inline Word addWords(Word &res, Word x, Word y, Word carry) { Word a = x + y; Word cy = a < x; Word b = a + carry; Word cz = b < a; res = b; return cy + cz; } inline Word addWords(Word &res, Word x, Word y) { Word a = x + y; res = a; return a < x; } Word add(DstPtr resp, SrcPtr xp, SrcPtr yp, Word xyn, Word carry = 0) { srep(i, xyn) carry = addWords(resp[i], xp[i], yp[i], carry); return carry; } Word add1(DstPtr resp, SrcPtr xp, Word xn, Word y) { srep(i, xn) { if(y == 0) return 0; Word a = xp[i] + y; resp[i] = a; y = a < y; } return y; } Word add(DstPtr resp, SrcPtr xp, Word xn, SrcPtr yp, Word yn) { if(xn < yn) swap(xp, yp), swap(xn, yn); Word c = add(resp, xp, yp, yn); return add1(resp + yn, xp + yn, xn - yn, c); } Word addMul(SrcDstPtr resp, SrcPtr xp, Size xn, Word y) { Word carry = 0; srep(i, xn) { Word lo, hi; multiplyWords(lo, hi, xp[i], y); carry = hi + addWords(resp[i], resp[i], lo, carry); } return carry; } void multiply(DstPtr resp, SrcPtr xp, Size xn, SrcPtr yp, Size yn) { assert(resp != xp && resp != yp); Size rn = xn + yn; srep(i, rn) resp[i] = 0; srep(i, yn) resp[i + xn] += addMul(resp + i, xp, xn, yp[i]); } typedef unsigned Num; struct ModNum { Num x; ModNum() : x() {} explicit ModNum(Num x_) : x(x_) {} Num get() const { return x; } }; typedef const ModNum *SrcPoly; typedef ModNum *DstPoly, SrcDstPoly; static const int NumBits = 32; template struct TmpBuffer { T *buf; TmpBuffer(Size n) { buf = new T[n]; } ~TmpBuffer() { delete[] buf; } T *data() { return buf; } }; typedef TmpBuffer TmpWords; typedef TmpBuffer TmpPoly; struct ModData { ModData() : mod() {} explicit ModData(Num m) { init(m); } Num getMod() const { return mod; } ModNum getOne() const { return ModNum(1); } void init(Num m) { if(NumBits != 32) cerr << "ModData: unimplemented" << endl, abort(); assert(m > 1); mod = m; shift = NumBits - getBitSize(m); den = m << shift; iden = reciprocalNum(den); } ModNum reduce(Num lo) const { return reduceSmall(lo, Num()); } //hi:lo < m:0 ModNum reduceSmall(Num lo, Num hi) const { Num quot, rem; div2by1ShiftSmall(hi, lo, quot, rem); return ModNum(rem); } ModNum mulMod(ModNum x, ModNum y) const { unsigned long long z = (unsigned long long)x.x * y.x; return reduceSmall((unsigned)z, (unsigned)(z >> 32)); } ModNum addMod(ModNum x, ModNum y) const { Num z = x.x + y.x; bool b = z >= mod || z < x.x; return ModNum(z - (mod & (-(Num)b))); } ModNum subtractMod(ModNum x, ModNum y) const { Num z = x.x + (mod - y.x); bool b = z >= mod || z < x.x; return ModNum(z - (mod & (-(Num)b))); } ModNum negateMod(ModNum x) const { return ModNum(x.x == 0 ? 0 : mod - x.x); } bool equal(ModNum x, ModNum y) const { return x.x == y.x; } bool equalModData(ModData that) const { return mod == that.mod; } int bitSize() const { return 32 - shift; } private: Num mod; int shift; Num den, iden; void div2by1(Num numHigh, Num numLow, Num ", Num &rem) const { Num qh, ql, r, mask; unsigned long long q = (unsigned long long)numHigh * iden + numLow; qh = (unsigned)(q >> 32), ql = (unsigned)q; qh += numHigh + 1; r = numLow - qh * den; mask = -(r > ql); qh += mask; r += den & mask; if(r >= den) { qh += 1; r -= den; } quot = qh; rem = r; } void div2by1ShiftSmall(Num numHigh, Num numLow, Num ", Num &rem) const { if(shift == 0) { div2by1(numHigh, numLow, quot, rem); } else { div2by1(numHigh << shift | numLow >> (NumBits - shift), numLow << shift, quot, rem); rem >>= shift; } } static Num reciprocalNum(Num d) { return (unsigned)(0xffffffffffffffffULL / d); } }; void packToBits(DstPtr resp, Size resn, SrcPoly xp, Size xn, int width, int stride) { if(width > WordBits) { cerr << "unimplemented" << endl; abort(); } srep(i, resn) resp[i] = 0; Size pos = 0; int sa = stride / WordBits, sb = stride % WordBits; int offset = 0; srep(i, xn) { assert(pos < resn); Word 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, ModData mod) { if(NumBits != 32) { cerr << "unimplemented" << endl; abort(); } assert(stride > 0); int K = (stride - 1) / NumBits + 1; TmpPoly bufBases(K); ModNum base = mod.reduceSmall(0, 1); ModNum *bases = bufBases.data(); bases[0] = mod.reduce(1); rep(i, K - 1) bases[i + 1] = mod.mulMod(bases[i], base); int pos = 0, offset = 0; srep(i, resn) { ModNum cur = ModNum(); for(int k = 0; k < K; ++ k) { int smallwidth = k != K - 1 ? NumBits : stride - NumBits * k; Word mask = ((Word)2 << (smallwidth - 1)) - 1; if(pos >= xn) { 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))); num &= mask; cur = mod.addMod(cur, mod.mulMod(mod.reduce(num), bases[k])); if(offset >= WordBits) offset -= WordBits, ++ pos; } resp[i] = cur; } } void multiplyPolynomialsKS(DstPoly resp, Size resn, SrcPoly xp, Size xn, SrcPoly yp, Size yn, ModData mod) { int width = mod.bitSize(); int stride = width * 2 + getBitSize(max(xn, yn)); Size tmpxn = static_cast((BitSize)xn * stride / WordBits + 1); Size tmpyn = static_cast((BitSize)yn * stride / WordBits + 1); Size tmpzn = tmpxn + tmpyn; TmpWords tmpx(tmpxn), tmpy(tmpyn), tmpz(tmpzn); packToBits(tmpx.data(), tmpxn, xp, xn, width, stride); packToBits(tmpy.data(), tmpyn, yp, yn, width, stride); multiply(tmpz.data(), tmpx.data(), tmpxn, tmpy.data(), tmpyn); unpackFromBitsMod(resp, resn, tmpz.data(), tmpzn, stride, mod); } void multiplyPolynomials(DstPoly resp, Size resn, SrcPoly xp, Size xn, SrcPoly yp, Size yn, ModData mod) { return multiplyPolynomialsKS(resp, resn, xp, xn, yp, yn, mod); } void squarePolynomial(DstPoly resp, Size resn, SrcPoly xp, Size xn, ModData mod) { return multiplyPolynomials(resp, resn, xp, xn, xp, xn, mod); } void addPolynomials(DstPoly resp, SrcPoly xp, SrcPoly yp, Size xyn, ModData mod) { srep(i, xyn) resp[i] = mod.addMod(xp[i], yp[i]); } void subtractPolynomials(DstPoly resp, SrcPoly xp, SrcPoly yp, Size xyn, ModData mod) { srep(i, xyn) resp[i] = mod.subtractMod(xp[i], yp[i]); } void multiplyConstantPolynomial(DstPoly resp, SrcPoly xp, Size xn, ModNum y, ModData mod) { srep(i, xn) resp[i] = mod.mulMod(xp[i], y); } void zeroPolynomial(DstPoly resp, Size resn) { srep(i, resn) resp[i] = ModNum(); } ModNum evaluatePolynomial(SrcPoly xp, Size xn, ModNum val, ModData mod) { if(xn == 0) return ModNum(); ModNum res = xp[xn - 1]; for(Size i = xn - 1; i > 0; -- i) res = mod.addMod(mod.mulMod(res, val), xp[i - 1]); return res; } void inversePowerSeries(DstPoly resp, Size resn, SrcPoly xp, Size xn, ModData mod) { if(resn == 0) return; assert(resp != xp); assert(mod.equal(xp[0], ModNum(1))); ModNum two = mod.reduce(2); TmpPoly tmp1(resn), tmp2(resn); zeroPolynomial(resp, resn); resp[0] = xp[0]; Size curn = 1; while(curn < resn) { Size nextn = min(resn, curn * 2); squarePolynomial(tmp1.data(), nextn, resp, curn, mod); multiplyPolynomials(tmp2.data(), nextn, tmp1.data(), nextn, xp, min(nextn, xn), mod); multiplyConstantPolynomial(resp, resp, curn, two, mod); subtractPolynomials(resp, resp, tmp2.data(), nextn, mod); curn = nextn; } } void reversePolynomial(DstPoly resp, SrcPoly xp, Size xn) { if(resp == xp) { reverse(resp, resp + xn); } else { srep(i, xn) resp[xn - 1 - i] = xp[i]; } } void precomputePolynomialInverse(DstPoly resp, Size resn, SrcPoly xp, Size xn, ModData mod) { TmpPoly tmp(xn); reversePolynomial(tmp.data(), xp, xn); inversePowerSeries(resp, resn, tmp.data(), xn, mod); } void dividePolynomialPrecomputed(DstPoly resp, Size resn, SrcPoly revxp, Size xn, SrcPoly invp, ModData mod) { //resn = xn - yn + 1 //size of inv >= resn multiplyPolynomials(resp, resn, revxp, xn, invp, resn, mod); reversePolynomial(resp, resp, resn); } void remainderPolynomialPrecomputed(DstPoly resp, SrcPoly xp, Size xn, SrcPoly yp, Size yn, SrcPoly invp, ModData mod) { if(xn < yn) { srep(i, xn) resp[i] = xp[i]; zeroPolynomial(resp + xn, yn - 1 - xn); return; } if(yn == 0) abort(); if(yp[yn - 1].get() != 1) { cerr << "unimplemented" << endl; abort(); } if(yn == 1) return; Size quotn = xn - yn + 1; TmpPoly revxp(xn), quot(quotn), quotmul(yn - 1); reversePolynomial(revxp.data(), xp, xn); dividePolynomialPrecomputed(quot.data(), quotn, revxp.data(), xn, invp, mod); multiplyPolynomials(quotmul.data(), yn - 1, quot.data(), min(quotn, yn - 1), yp, yn - 1, mod); subtractPolynomials(resp, xp, quotmul.data(), yn - 1, mod); } } //namespace multiprecision namespace multiprecisionxx { using namespace multiprecision; struct Polynomial { ModData mod; vector coeffs; Polynomial() {} Polynomial(ModData mod_) { initMod(mod_); } Polynomial(ModData mod_, Size initsize): mod(mod_), coeffs(initsize) { } Size size() const { return (Size)coeffs.size(); } ModNum *data() { return coeffs.empty() ? NULL : coeffs.data(); } const ModNum *data() const { return coeffs.empty() ? NULL : coeffs.data(); } Polynomial getOne() const { Polynomial res(mod, 1); res.setCoeff(0, mod.getOne()); return res; } void initMod(ModData mod_) { mod = mod_; coeffs.clear(); } void setCoeff(Size i, ModNum x) { if(i >= size()) resize(i+1); coeffs[i] = x; normalize(); } ModNum getCoeff(Size i) const { if(i < 0 || i >= size()) return ModNum(); else return coeffs[i]; } Polynomial operator*(const Polynomial &that) const { Polynomial res(mod, size() + that.size() - 1); if(this == &that) { squarePolynomial(res.data(), res.size(), data(), size(), mod); } else { assert(compatibleTo(that)); multiplyPolynomials(res.data(), res.size(), data(), size(), that.data(), that.size(), mod); } res.normalize(); return res; } Polynomial precomputeInverse(Size n) const { Polynomial res(mod, n); precomputePolynomialInverse(res.data(), res.size(), data(), size(), mod); return res; } Polynomial remainder(const Polynomial &that, const Polynomial &inv) const { assert(compatibleTo(that)); assert(size() < that.size() || size() - that.size() + 1 <= inv.size()); Polynomial res(mod, that.size() - 1); remainderPolynomialPrecomputed(res.data(), data(), size(), that.data(), that.size(), inv.data(), mod); res.normalize(); return res; } Polynomial powMod(long long k, const Polynomial &mod, const Polynomial &inv) const { assert(compatibleTo(mod) && compatibleTo(inv)); assert(k >= 0); if(k == 0) return getOne(); return powModRec(*this, k, mod, inv); } void resize(Size n) { coeffs.resize(n); } void normalize() { while(!coeffs.empty() && coeffs.back().x == 0) coeffs.pop_back(); } bool compatibleTo(const Polynomial &that) const { return mod.equalModData(that.mod); } private: static Polynomial powModRec(const Polynomial &base, long long k, const Polynomial &mod, const Polynomial &inv) { assert(k >= 1); if(k == 1) return base; Polynomial p = powModRec(base, k / 2, mod, inv); p = p * p; p = p.remainder(mod, inv); if(k % 2 != 0) { p = p * base; p = p.remainder(mod, inv); } return p; } }; } namespace testtest { #pragma region tests struct Xor128 { unsigned x, y, z, w; Xor128() : x(123456789), y(362436069), z(521288629), w(88675123) {} //[0, 2^32) unsigned operator()() { unsigned t = x ^ (x << 11); x = y; y = z; z = w; return w = w ^ (w >> 19) ^ (t ^ (t >> 8)); } //[0, 2^64) unsigned long long nextLL() { unsigned x = (*this)(); unsigned y = (*this)(); return (unsigned long long)x << 32 | y; } //[0, n) unsigned operator()(unsigned n) { unsigned mask = calculateMask(n - 1), x; do { x = (*this)() & mask; } while(x >= n); return x; } //[0, n) signed int operator()(signed int n) { return (*this)((unsigned int)n); } //[L, U] signed int operator()(signed int L, signed int U) { return L + (*this)(U - L + 1); } //[0, n) unsigned long long operator()(unsigned long long n) { unsigned long long mask = calculateMask(n - 1), x; do { x = (*this).nextLL() & mask; } while(x >= n); return x; } //[0, n) signed long long operator()(signed long long n) { return (*this)((unsigned long long)n); } //[L, U] signed long long operator()(signed long long L, signed long long U) { return L + (*this)(U - L + 1); } private: static unsigned calculateMask(unsigned v) { v |= v >> 1; v |= v >> 2; v |= v >> 4; v |= v >> 8; v |= v >> 16; return v; } static unsigned long long calculateMask(unsigned long long v) { v |= v >> 1; v |= v >> 2; v |= v >> 4; v |= v >> 8; v |= v >> 16; v |= v >> 32; return v; } }; struct GModInt { static int Mod; unsigned x; GModInt() : x(0) {} GModInt(signed sig) { int sigt = sig % Mod; if(sigt < 0) sigt += Mod; x = sigt; } GModInt(signed long long sig) { int sigt = sig % Mod; if(sigt < 0) sigt += Mod; x = sigt; } int get() const { return (int)x; } GModInt &operator+=(GModInt that) { if((x += that.x) >= (unsigned)Mod) x -= Mod; return *this; } GModInt &operator-=(GModInt that) { if((x += Mod - that.x) >= (unsigned)Mod) x -= Mod; return *this; } GModInt &operator*=(GModInt that) { x = (unsigned long long)x * that.x % Mod; return *this; } GModInt &operator/=(GModInt that) { return *this *= that.inverse(); } GModInt operator+(GModInt that) const { return GModInt(*this) += that; } GModInt operator-(GModInt that) const { return GModInt(*this) -= that; } GModInt operator*(GModInt that) const { return GModInt(*this) *= that; } GModInt operator/(GModInt that) const { return GModInt(*this) /= that; } //Modと素であることが保証されるかどうか確認すること! GModInt inverse() const { signed a = x, b = Mod, u = 1, v = 0; while(b) { signed t = a / b; a -= t * b; std::swap(a, b); u -= t * v; std::swap(u, v); } if(u < 0) u += Mod; GModInt res; res.x = (unsigned)u; return res; } bool operator==(GModInt that) const { return x == that.x; } bool operator!=(GModInt that) const { return x != that.x; } GModInt operator-() const { GModInt t; t.x = x == 0 ? 0 : Mod - x; return t; } }; int GModInt::Mod = 0; struct Polynomial { typedef GModInt Coef; typedef Coef Val; vector coef; //... + coef[2] x^2 + coef[1] x + coef[0] Polynomial() {} explicit Polynomial(int n) : coef(n) {} static Polynomial One() { Polynomial r(1); r.coef[0] = 1; return r; } bool iszero() const { return coef.empty(); } int degree1() const { return (int)coef.size(); } //degree + 1 int resize(int d) { if(degree1() < d) coef.resize(d); return d; } const Coef operator[](int i) const { return i >= degree1() ? Coef() : coef[i]; } void canonicalize() { int i = (int)coef.size(); while(i > 0 && coef[i - 1] == Coef()) i --; coef.resize(i); } Val evalute(Val x) const { int d = degree1(); Val t = 0, y = 1; rep(i, d) { t += y * coef[i]; y *= x; } return t; } Polynomial &operator+=(const Polynomial &that) { int d = resize(that.degree1()); for(int i = 0; i < d; i ++) coef[i] += that[i]; canonicalize(); return *this; } Polynomial operator+(const Polynomial &that) const { return Polynomial(*this) += that; } Polynomial &operator-=(const Polynomial &that) { int d = resize(that.degree1()); for(int i = 0; i < d; i ++) coef[i] -= that[i]; canonicalize(); return *this; } Polynomial operator-(const Polynomial &that) const { return Polynomial(*this) -= that; } Polynomial operator-() const { int d = degree1(); Polynomial res(d); for(int i = 0; i < d; i ++) res.coef[i] = - coef[i]; return res; } //naive Polynomial operator*(const Polynomial &that) const { if(iszero() || that.iszero()) return Polynomial(); int x = degree1(), y = that.degree1(), d = x + y - 1; Polynomial res(d); rep(i, x) rep(j, y) res.coef[i + j] += coef[i] * that.coef[j]; res.canonicalize(); return res; } //long division pair divmod(const Polynomial &that) const { int x = degree1() - 1, y = that.degree1() - 1; int d = max(0, x - y); Polynomial q(d + 1), r = *this; for(int i = x; i >= y; i --) { Coef t = r.coef[i] / that.coef[y]; q.coef[i - y] = t; assert(t * that.coef[y] == r.coef[i]); r.coef[i] = 0; if(t == 0) continue; for(int j = 0; j < y; j ++) r.coef[i - y + j] -= t * that.coef[j]; } q.canonicalize(); r.canonicalize(); return make_pair(q, r); } Polynomial operator/(const Polynomial &that) const { return divmod(that).first; } Polynomial operator%(const Polynomial &that) const { return divmod(that).second; } }; namespace test_multiprecision { using namespace multiprecision; Xor128 xor128; void randomWords(DstPtr xp, Size xn) { srep(i, xn) { #ifdef MP_64 xp[i] = xor128.nextLL(); #else xp[i] = xor128(); #endif } } void randomPolynomial(DstPoly xp, Size xn, ModData mod) { srep(i, xn) xp[i] = ModNum(xor128(0, mod.getMod() - 1)); } const unsigned TestRandomPrime = 1000000007; unsigned testReduceRandomMod(SrcPtr xp, Size n) { const unsigned P = TestRandomPrime; #ifdef MP_64 unsigned x = ((unsigned long long)((1ULL << 32) % P) << 32) % P; #else unsigned x = (1ULL << 32) % P; #endif unsigned r = 0; for(Size i = n; i > 0; -- i) r = ((unsigned long long)r * x + xp[i - 1] % P) % P; return r; } const int NumTests = 100000; Num testRandomMod(int maxbits = 32) { if(NumBits != 32) abort(); int l = xor128(3, maxbits); return xor128(2, (2U << (l - 1)) - 1); } void testMod() { if(NumBits != 32) abort(); rep(ii, NumTests) { unsigned m = testRandomMod(); ModData mod(m); unsigned x = xor128(0, m - 1); unsigned y = xor128(0, m - 1); unsigned long long z = (unsigned long long)x * y; unsigned a = mod.reduceSmall((unsigned)z, (unsigned)(z >> 32)).get(); unsigned b = z % m; assert(a == b); unsigned c = mod.mulMod(ModNum(x), ModNum(y)).get(); assert(c == b); unsigned d = mod.addMod(ModNum(x), ModNum(y)).get(); assert(d == ((unsigned long long)x + y) % m); unsigned e = mod.subtractMod(ModNum(x), ModNum(y)).get(); assert(e == ((unsigned long long)x + m - y) % m); } } void testMultiply() { const int MaxN = 20; TmpWords xs(MaxN), ys(MaxN), res(MaxN * 2); rep(ii, NumTests) { int xn = xor128(1, MaxN), yn = xor128(1, MaxN); randomWords(xs.data(), xn); randomWords(ys.data(), yn); unsigned xa = testReduceRandomMod(xs.data(), xn); unsigned ya = testReduceRandomMod(ys.data(), yn); multiply(res.data(), xs.data(), xn, ys.data(), yn); unsigned resa = testReduceRandomMod(res.data(), xn + yn); unsigned ansa = (unsigned long long)xa * ya % TestRandomPrime; assert(ansa == resa); } } void testMultiplyKS() { const int MaxN = 20; TmpPoly xs(MaxN), ys(MaxN), res(MaxN * 2); rep(ii, NumTests) { Num m = testRandomMod(); ModData mod(m); int xn = xor128(1, MaxN), yn = xor128(1, MaxN); randomPolynomial(xs.data(), xn, mod); randomPolynomial(ys.data(), yn, mod); ModNum randomVal = ModNum(xor128(0, m - 1)); ModNum xa = evaluatePolynomial(xs.data(), xn, randomVal, mod); ModNum ya = evaluatePolynomial(ys.data(), yn, randomVal, mod); res.data()[xn + yn - 1].x = 0x12345678; multiplyPolynomialsKS(res.data(), xn + yn - 1, xs.data(), xn, ys.data(), yn, mod); assert(res.data()[xn + yn - 1].x == 0x12345678); ModNum resa = evaluatePolynomial(res.data(), xn + yn - 1, randomVal, mod); ModNum ansa = mod.mulMod(xa, ya); assert(ansa.x == resa.x); } } void testRemainderPolynomial() { const int MaxN = 20; TmpPoly xs(MaxN), ys(MaxN), invy(MaxN), res(MaxN); rep(ii, NumTests) { Num m = testRandomMod(31); ModData mod(m); int xn = xor128(1, MaxN), yn = xor128(1, xn); randomPolynomial(xs.data(), xn, mod); randomPolynomial(ys.data(), yn - 1, mod); ys.data()[yn - 1] = mod.reduce(1); precomputePolynomialInverse(invy.data(), xn - yn + 1, ys.data(), yn, mod); remainderPolynomialPrecomputed(res.data(), xs.data(), xn, ys.data(), yn, invy.data(), mod); GModInt::Mod = m; Polynomial polyx(xn), polyy(yn); rep(i, xn) polyx.coef[i].x = xs.data()[i].x; rep(i, yn) polyy.coef[i].x = ys.data()[i].x; Polynomial polyres = polyx % polyy; rep(i, yn - 1) assert(polyres[i].x == res.data()[i].x); } } } void tests() { using namespace test_multiprecision; testMod(); testMultiply(); testMultiplyKS(); testRemainderPolynomial(); cerr << "tests: completed" << endl; } #pragma endregion } 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(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; } }; typedef ModInt<1000000007> mint; vector doDP(const int a[6], int n) { int maxX = a[5] * n; vector > dp(n+1, vector(maxX+1)); dp[0][0] = 1; rep(k, 6) { int t = a[k]; for(int i = n; i >= 0; -- i) { int maxx = k == 0 ? 0 : i * a[k-1]; rer(j, 0, maxx) { mint x = dp[i][j]; if(x.get() == 0) continue; rer(l, 1, n-i) dp[i+l][j + l * t] += 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 multiprecisionxx; int n = (int)c.size(); ModData mod = ModData(mint::Mod); Polynomial poly(mod); poly.setCoeff(n, mod.getOne()); rep(i, n) poly.setCoeff(i, mod.negateMod(ModNum(c[i].x))); Polynomial inv = poly.precomputeInverse(n + 1); Polynomial x(mod); x.setCoeff(1, mod.getOne()); Polynomial p = x.powMod(N, poly, inv); vector res(n); rep(i, n) res[i].x = p.getCoeff(i).get(); return res; } int main() { // testtest::tests(); long long N; int P; int C; while(~scanf("%lld%d%d", &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 = doDP(primes, P); vector cntC = doDP(composites, C); int X = P * primes[5] + C * composites[5]; vector ways(X+1); rep(i, cntP.size()) rep(j, cntC.size()) ways[i + j] += cntP[i] * cntC[j]; vector c(X); rep(i, X) c[i] = ways[X - i]; vector coeffs = computeLinearRecurrentSequenceCoefficients(c, X + (N-1)); mint ans = accumulate(all(coeffs), mint()); printf("%d\n", ans.get()); } return 0; }