#51474 (C++11) No.214 素数サイコロと合成数サイコロ (3-Medium)

提出ソース
結果

問題	No.214 素数サイコロと合成数サイコロ (3-Medium)
ユーザー	anta
提出日時	2015-09-29 01:40:50
言語	C++11 (gcc 11.4.0)
結果	AC
実行時間	759 ms / 3,000 ms
コード長	26,034 bytes
コンパイル時間	1,854 ms
コンパイル使用メモリ	107,848 KB
実行使用メモリ	4,380 KB
最終ジャッジ日時	2023-09-26 17:29:20
合計ジャッジ時間	5,244 ms
ジャッジサーバーID （参考情報）	judge11 / judge15
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	743 ms 4,380 KB
testcase_01	AC	702 ms 4,380 KB
testcase_02	AC	759 ms 4,376 KB
権限があれば一括ダウンロードができます
ソースコード

raw source code
#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 &quot, 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 &quot, 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;
}
yukicoder

結果

テストケース

ソースコード
