結果
| 問題 |
No.214 素数サイコロと合成数サイコロ (3-Medium)
|
| コンテスト | |
| ユーザー |
anta
|
| 提出日時 | 2015-09-29 01:40:50 |
| 言語 | C++11(廃止可能性あり) (gcc 13.3.0) |
| 結果 |
AC
|
| 実行時間 | 819 ms / 3,000 ms |
| コード長 | 26,034 bytes |
| コンパイル時間 | 1,449 ms |
| コンパイル使用メモリ | 111,576 KB |
| 実行使用メモリ | 5,376 KB |
| 最終ジャッジ日時 | 2024-07-19 11:28:02 |
| 合計ジャッジ時間 | 5,144 ms |
|
ジャッジサーバーID (参考情報) |
judge3 / judge2 |
(要ログイン)
| ファイルパターン | 結果 |
|---|---|
| other | AC * 3 |
ソースコード
#include <string>
#include <vector>
#include <algorithm>
#include <numeric>
#include <set>
#include <map>
#include <queue>
#include <iostream>
#include <sstream>
#include <cstdio>
#include <cmath>
#include <ctime>
#include <cstring>
#include <cctype>
#include <cassert>
#include <limits>
#include <functional>
#ifndef __GNUC__
#include <intrin.h>
#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<int> vi; typedef pair<int,int> pii; typedef vector<pair<int,int> > vpii; typedef long long ll;
template<typename T, typename U> inline void amin(T &x, U y) { if(y < x) x = y; }
template<typename T, typename U> 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<Word>(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<int>(res) + 1;
#else
return 64 - __builtin_clzll(x);
#endif
#else
#ifndef __GNUC__
unsigned long res;
_BitScanReverse(&res, x);
return static_cast<int>(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<typename T>
struct TmpBuffer {
T *buf;
TmpBuffer(Size n) {
buf = new T[n];
}
~TmpBuffer() {
delete[] buf;
}
T *data() { return buf; }
};
typedef TmpBuffer<Word> TmpWords;
typedef TmpBuffer<ModNum> 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<Word>(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<Num>(xp[pos] >> offset);
offset += smallwidth;
if(pos + 1 < xn && offset > WordBits)
num |= static_cast<Num>(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<Size>((BitSize)xn * stride / WordBits + 1);
Size tmpyn = static_cast<Size>((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<ModNum> 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; //... + 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<Polynomial, Polynomial> 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<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(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<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] = 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<mint> computeLinearRecurrentSequenceCoefficients(const vector<mint> &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<mint> 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<mint> cntP = doDP(primes, P);
vector<mint> cntC = doDP(composites, C);
int X = P * primes[5] + C * composites[5];
vector<mint> ways(X+1);
rep(i, cntP.size()) rep(j, cntC.size())
ways[i + j] += cntP[i] * cntC[j];
vector<mint> c(X);
rep(i, X)
c[i] = ways[X - i];
vector<mint> coeffs = computeLinearRecurrentSequenceCoefficients(c, X + (N-1));
mint ans = accumulate(all(coeffs), mint());
printf("%d\n", ans.get());
}
return 0;
}
anta