結果

問題 No.215 素数サイコロと合成数サイコロ (3-Hard)
コンテスト
ユーザー anta
提出日時 2015-10-06 04:34:44
言語 C++11（廃止可能性あり）
(gcc 13.3.0)
結果
CE  
(最新)
AC  
(最初)
実行時間 -
コード長 38,254 bytes
コンパイル時間 710 ms
コンパイル使用メモリ 84,776 KB
最終ジャッジ日時 2024-11-14 19:18:16
合計ジャッジ時間 1,168 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge1
このコードへのチャレンジ
（要ログイン）
コンパイルエラー時のメッセージ・ソースコードは、提出者また管理者しか表示できないようにしております。（リジャッジ後のコンパイルエラーは公開されます）
ただし、clay言語の場合は開発者のデバッグのため、公開されます。

コンパイルメッセージ
main.cpp:562:8: error: ‘__m128i’ does not name a type
  562 | inline __m128i add_mod_128i(const __m128i &x, const __m128i &y, const __m128i &mod) {
      |        ^~~~~~~
main.cpp:570:8: error: ‘__m128i’ does not name a type
  570 | inline __m128i sub_mod_128i(const __m128i &x, const __m128i &y, const __m128i &zero, const __m128i &mod) {
      |        ^~~~~~~
main.cpp:768:1: warning: inline variables are only available with ‘-std=c++17’ or ‘-std=gnu++17’
  768 | inline unsigned getpart_epi32(__m128i x, const int i) {
      | ^~~~~~
main.cpp:768:31: error: ‘__m128i’ was not declared in this scope
  768 | inline unsigned getpart_epi32(__m128i x, const int i) {
      |                               ^~~~~~~
main.cpp:768:42: error: expected primary-expression before ‘const’
  768 | inline unsigned getpart_epi32(__m128i x, const int i) {
      |                                          ^~~~~
main.cpp:768:53: error: expression list treated as compound expression in initializer [-fpermissive]
  768 | inline unsigned getpart_epi32(__m128i x, const int i) {
      |                                                     ^

ソースコード

diff # 

#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>
#include <cstdint>

#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

#if defined(_M_X64) || defined(__amd64__)
#define MP_64
#endif


struct MontgomeryModInt {
	static const int Mod = 1000000007;
	static const int W = 32;
	static const int R = (1ULL << W) % Mod;	//= 2^W
	static const int InvR = 518424770;		//= R^{-1} (mod Mod)
	static const int InvNegMod = (int)2226617417U;	//= (-Mod)^{-1} (mod R)
#if defined(_MSC_VER) || __cplusplus > 199711L
	static_assert(InvR < Mod && ((unsigned long long)R * InvR) % Mod == 1, "InvR = R^{-1} (mod Mod)");
	static_assert((((unsigned long long)((1ULL << W) - Mod) * (unsigned long long)(unsigned)InvNegMod) & 0xffffffffULL) == 1, "InvNegMod = (-Mod)^{-1} (mod R)");
#endif

	unsigned residue;	//residue = value * R % Mod

	MontgomeryModInt(): residue(0) { }
	MontgomeryModInt(signed sig): residue(multR(sig)) { }
	MontgomeryModInt(signed long long sig): residue(multR((signed)(sig % Mod))) { }

	static unsigned multR(signed sig) {
		signed sigt = (signed long long)sig * R % Mod;
		if(sigt < 0) sigt += Mod;
		return (unsigned)sigt;
	}

	int get() const { return (int)reduce(residue); }

	static unsigned reduce(unsigned T) {
		unsigned m = T * (unsigned)InvNegMod;
		unsigned t = (unsigned)((T + (unsigned long long)m * Mod) >> W);
		if(t >= Mod) t -= Mod;
		return t;
	}
	static unsigned reduce(unsigned long long T) {
		unsigned m = (unsigned)T * (unsigned)InvNegMod;
		unsigned t = (unsigned)((T + (unsigned long long)m * Mod) >> W);
		if(t >= Mod) t -= Mod;
		return t;
	}

	MontgomeryModInt &operator+=(MontgomeryModInt that) {
		if((residue += that.residue) >= Mod) residue -= Mod;
		return *this;
	}
	MontgomeryModInt &operator-=(MontgomeryModInt that) {
		if((residue += Mod - that.residue) >= Mod) residue -= Mod;
		return *this;
	}
	MontgomeryModInt &operator*=(MontgomeryModInt that) {
		residue = reduce((unsigned long long)residue * that.residue);
		return *this;
	}
	MontgomeryModInt &operator/=(MontgomeryModInt that) {
		return *this *= that.inverse();
	}

	MontgomeryModInt operator+(MontgomeryModInt that) const { return MontgomeryModInt(*this) += that; }
	MontgomeryModInt operator-(MontgomeryModInt that) const { return MontgomeryModInt(*this) -= that; }
	MontgomeryModInt operator*(MontgomeryModInt that) const {
		MontgomeryModInt res;
		res.residue = reduce((unsigned long long)residue * that.residue);
		return res;
	}
	MontgomeryModInt operator/(MontgomeryModInt that) const { return MontgomeryModInt(*this) /= that; }

	MontgomeryModInt inverse() const {
		//a^{-1} R = (a R)^{-1} * R^2
		static const int SquaredR = (unsigned long long)R * R % Mod;
		signed a = residue, 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;
		MontgomeryModInt res;
		res.residue = (unsigned long long)(unsigned)u * SquaredR % Mod;
		return res;
	}

	MontgomeryModInt operator-() const { MontgomeryModInt t; t.residue = residue == 0 ? 0 : Mod - residue; return t; }

	bool operator==(const MontgomeryModInt &that) const { return residue == that.residue; }
	bool operator!=(const MontgomeryModInt &that) const { return residue != that.residue; }
};
typedef MontgomeryModInt mint;

namespace mod_polynomial {

struct Polynomial {
	typedef mint R;
	static R ZeroR() { return R(); }
	static R OneR() { return R(1); }
	static bool IsZeroR(R r) { return r == ZeroR(); }

	std::vector<R> coefs;

	Polynomial() { }
	explicit Polynomial(R c0): coefs(1, c0) { }
	explicit Polynomial(R c0, R c1): coefs(2) { coefs[0] = c0, coefs[1] = c1; }
	template<typename It> Polynomial(It be, It en): coefs(be, en) { }

	static Polynomial Zero() { return Polynomial(); }
	static Polynomial One() { return Polynomial(OneR()); }
	static Polynomial X() { return Polynomial(ZeroR(), OneR()); }

	void resize(int n) { coefs.resize(n); }
	void clear() { coefs.clear(); }

	R *data() { return coefs.empty() ? NULL : &coefs[0]; }
	const R *data() const { return coefs.empty() ? NULL : &coefs[0]; }

	int size() const { return static_cast<int>(coefs.size()); }
	bool empty() const { return coefs.empty(); }
	int degree() const { return size() - 1; }

	bool normalized() const { return coefs.empty() || coefs.back() != ZeroR(); }
	bool monic() const { return !coefs.empty() && coefs.back() == OneR(); }

	R get(int i) const { return 0 <= i && i < size() ? coefs[i] : ZeroR(); }

	void set(int i, R x) {
		if(size() <= i)
			resize(i + 1);
		coefs[i] = x;
	}

	void normalize() { while(!empty() && IsZeroR(coefs.back())) coefs.pop_back(); }

	R evaluate(R x) const {
		if(empty()) return R();
		R r = coefs.back();
		for(int i = size() - 2; i >= 0; -- i) {
			r *= x;
			r += coefs[i];
		}
		return r;
	}

	Polynomial &operator+=(const Polynomial &that) {
		int m = size(), n = that.size();
		if(m < n) resize(n);
		_add(data(), that.data(), n);
		return *this;
	}
	Polynomial operator+(const Polynomial &that) const {
		return Polynomial(*this) += that;
	}

	Polynomial &operator-=(const Polynomial &that) {
		int m = size(), n = that.size();
		if(m < n) resize(n);
		_subtract(data(), that.data(), n);
		return *this;
	}
	Polynomial operator-(const Polynomial &that) const {
		return Polynomial(*this) -= that;
	}

	Polynomial &operator*=(R r) {
		_multiply_1(data(), size(), r);
		return *this;
	}
	Polynomial operator*(R r) const {
		Polynomial res;
		res.resize(size());
		_multiply_1(res.data(), data(), size(), r);
		return res;
	}
	Polynomial operator*(int z) const {
		if(z < 0) return -(*this * -z);
		if(z == 0) return Polynomial();
		if(z == 1) return *this;
		Polynomial res;
		res.resize(size());
		_multiply_1_int(res.data(), data(), size(), z);
		return res;
	}
	Polynomial operator/(int z) const {
		if(z < 0) return -(*this / -z);
		if(z == 0) abort();
		if(z == 1) return *this;
		Polynomial res;
		res.resize(size());
		_divide_1_int(res.data(), data(), size(), z);
		return res;
	}

	Polynomial operator*(const Polynomial &that) const {
		Polynomial r;
		multiply(r, *this, that);
		return r;
	}

	Polynomial &operator*=(const Polynomial &that) {
		multiply(*this, *this, that);
		return *this;
	}

	static void multiply(Polynomial &res, const Polynomial &p, const Polynomial &q) {
		int pn = p.size(), qn = q.size();

		if(pn < qn)
			return multiply(res, q, p);

		if(&res == &p || &res == &q) {
			Polynomial tmp;
			multiply(tmp, p, q);
			res = tmp;
			return;
		}

		if(qn == 0) {
			res.coefs.clear();
		}else {
			res.resize(pn + qn - 1);
			_multiply_select_method(res.data(), p.data(), pn, q.data(), qn);
		}
	}

	Polynomial operator<<(int shift) const {
		if(shift < 0) return *this >> (-shift);
		Polynomial res;
		res.resize(size() + shift);
		_copy(res.data() + shift, data(), size());
		return res;
	}

	Polynomial operator>>(int shift) const {
		if(shift < 0) return *this << (-shift);
		if(shift >= size()) return Polynomial();
		Polynomial res;
		res.resize(size() - shift);
		_copy(res.data(), data() + shift, size() - shift);
		return res;
	}

	Polynomial operator-() const {
		Polynomial res;
		res.resize(size());
		_negate(res.data(), data(), size());
		return res;
	}

	Polynomial precomputeInverse(int n) const {
		Polynomial res;
		res.resize(n);
		_precompute_inverse(res.data(), n, data(), size());
		return res;
	}

	static void divideRemainderPrecomputedInverse(Polynomial &quot, Polynomial &rem, const Polynomial &p, const Polynomial &q, const Polynomial &inv) {
		assert(&quot != &p && &quot != &q && &quot != &inv);
		int pn = p.size(), qn = q.size();
		assert(inv.size() >= pn - qn + 1);
		quot.resize(std::max(0, pn - qn + 1));
		rem.resize(qn - 1);
		_divide_remainder_precomputed_inverse(quot.data(), rem.data(), p.data(), pn, q.data(), qn, inv.data());
		quot.normalize();
		rem.normalize();
	}

	Polynomial computeRemainderPrecomputedInverse(const Polynomial &q, const Polynomial &inv) const {
		Polynomial quot, rem;
		divideRemainderPrecomputedInverse(quot, rem, *this, q, inv);
		return rem;
	}

	Polynomial powerMod(long long K, const Polynomial &q) const {
		int qn = q.size();
		assert(K >= 0 && qn > 0);
		assert(q.monic());
		if(qn == 1) return Polynomial();
		if(K == 0) return One();
		Polynomial inv = q.precomputeInverse(std::max(size() - qn + 1, qn));
		Polynomial p = this->computeRemainderPrecomputedInverse(q, inv);
		int l = 0;
		while((K >> l) > 1) ++ l;
		Polynomial res = p;
		for(-- l; l >= 0; -- l) {
			res *= res;
			res = res.computeRemainderPrecomputedInverse(q, inv);
			if(K >> l & 1) {
				res *= p;
				res = res.computeRemainderPrecomputedInverse(q, inv);
			}
		}
		return res;
	}

	static int TOOM22_THRESHOLD;
	static int TOOM33_THRESHOLD;
	static int TOOM44_THRESHOLD;
	static int SCHONHAGE_STRASSEN_THRESHOLD;

//private:
	class WorkSpaceStack;

	static void _fill_zero(R *res, int n);
	static void _copy(R *res, const R *p, int n);
	static void _negate(R *res, const R *p, int n);

	static void _add(R *p, const R *q, int n);
	static void _add(R *res, const R *p, int pn, const R *q, int qn);
	static void _subtract(R *p, const R *q, int n);
	static void _subtract(R *res, const R *p, int pn, const R *q, int qn);

	static void _multiply_select_method(R *res, const R *p, int pn, const R *q, int qn);
	static void _square_select_method(R *res, const R *p, int pn);
public:
	static int _schonhage_strassen_determine_first_k(int n);
	static int _schonhage_strassen_determine_recurse_k(int n, int k);
//private:

	static void _multiply_1(R *p, const R *q, int n, R c0);
	static void _multiply_1(R *p, int n, R c0);
	static void _add_multiply_1(R *p, const R *q, int n, R c0);

	static void _multiply_1_int(R *res, const R *p, int n, int z);
	static void _divide_1_int(R *res, const R *p, int n, int z);
	static void _multiply_power_of_two(R *res, const R *p, int n, int k);
	static void _divide_power_of_two(R *res, const R *p, int n, int k);

	static void _schoolbook_multiplication(R *res, const R *p, int pn, const R *q, int qn);

	static void _toom22(R *res, const R *p, int pn, const R *q, int qn, WorkSpaceStack &wss);
	static int _toom22_estimate_workspace(int pn, int qn);
	static void _toom33(R *res, const R *p, int pn, const R *q, int qn);
	static void _toom44(R *res, const R *p, int pn, const R *q, int qn);

	static void _schonhage_strassen(R *res, const R *p, int pn, const R *q, int qn);
	static void _schonhage_strassen_recurse(R *res, const R *X, const R *Y, int N, int k);
	static void _schonhage_strassen_decompose(R *x, const R *X, int k, int n, int m, R *tmp);
	static void _fft_make_omega_table(int n, int k, unsigned short *omega_table);
	static void _fft_rec(R *x, int n, int inc, int k, R *tmp, const unsigned short *omega_table);
	static void _fft_inverse_rec(R *x, int n, int k, R *tmp);

	static void _negacyclic_shift(R *res, const R *x, int n, int shift);
	static void _negacyclic_shift_add(R *p, int pn, const R *q, int qn, int shift);

	static void _reverse(R *res, const R *p, int pn);
	static void _inverse_power_series(R *res, int resn, const R *p, int pn);
	static void _precompute_inverse(R *res, int resn, const R *p, int pn);
	static void _divide_precomputed_inverse(R *res, int resn, const R *revp, int pn, const R *inv);
	static void _divide_remainder_precomputed_inverse(R *quot, R *rem, const R *p, int pn, const R *q, int qn, const R *inv);

	class WorkSpaceStack {
	public:
		explicit WorkSpaceStack(int n_): n(n_), used(0), max_used(0) {
			p = static_cast<R*>(std::malloc(n * sizeof(R)));
		}
		~WorkSpaceStack() {
			std::free(p);
		}

		R *allocate(int size) {
			if(size > n - used) {
				std::abort();
			}
			R *t = p + used;
			used += size;
			if(used > max_used) max_used = used;
			return t;
		}
		void deallocate(R *t, int size) {
			if(size > used || p + (used - size) != t) {
				std::abort();
			}
			used -= size;
		}

	private:
		R *p; int n, used, max_used;

		WorkSpaceStack(const WorkSpaceStack &);
		WorkSpaceStack &operator=(const WorkSpaceStack &);
	};

	class WorkSpace {
	public:
		explicit WorkSpace(WorkSpaceStack &stk_, int n_, bool fill_zero = false): stk(&stk_), n(n_) {
			p = stk->allocate(n_);
			if(fill_zero) _fill_zero(p, n);
		}
		explicit WorkSpace(int n_, bool fill_zero = false): stk(NULL), n(n_) {
			p = static_cast<R*>(std::malloc(n * sizeof(R)));
			if(fill_zero) _fill_zero(p, n);
		}
		~WorkSpace() {
			if(stk != NULL) {
				stk->deallocate(p, n);
			}else {
				std::free(p);
			}
		}
		R *get() { return p; }
	private:
		WorkSpaceStack *stk;
		R *p; int n;

		WorkSpace(const WorkSpace &);
		WorkSpace &operator=(const WorkSpace &);
	};
};
int Polynomial::TOOM22_THRESHOLD = 16;
int Polynomial::TOOM33_THRESHOLD = 120;
int Polynomial::TOOM44_THRESHOLD = 120;
int Polynomial::SCHONHAGE_STRASSEN_THRESHOLD = 120;

void Polynomial::_fill_zero(R *res, int n) {
	for(int i = 0; i < n; ++ i)
		res[i] = ZeroR();
}

void Polynomial::_copy(R *res, const R *p, int n) {
	for(int i = 0; i < n; ++ i)
		res[i] = p[i];
}

#ifndef __GNUC__

inline __m128i add_mod_128i(const __m128i &x, const __m128i &y, const __m128i &mod) {
	__m128i sum = _mm_add_epi32(x, y);
	__m128i lt_mod_mask = _mm_cmpgt_epi32(mod, sum);
	__m128i mod_sub = _mm_andnot_si128(lt_mod_mask, mod);
	__m128i result = _mm_sub_epi32(sum, mod_sub);
	return result;
}

inline __m128i sub_mod_128i(const __m128i &x, const __m128i &y, const __m128i &zero, const __m128i &mod) {
	__m128i diff = _mm_sub_epi32(x, y);
	__m128i lt_0_mask = _mm_cmpgt_epi32(zero, diff);
	__m128i mod_add = _mm_and_si128(lt_0_mask, mod);
	__m128i result = _mm_add_epi32(diff, mod_add);
	return result;
}

void Polynomial::_negate(R *res, const R *p, int n) {
	if(n >= 4) {
		__m128i zero = _mm_setzero_si128();
		__m128i mod = _mm_set1_epi32(R::Mod);
		for(int i = 0; i + 3 < n; i += 4) {
			__m128i x = _mm_loadu_si128(reinterpret_cast<const __m128i*>(p + i));
			__m128i result = sub_mod_128i(zero, x, zero, mod);
			_mm_storeu_si128(reinterpret_cast<__m128i*>(res + i), result);
		}
	}
	for(int i = n - (n & 3); i < n; ++ i)
		res[i] = -p[i];
}

void Polynomial::_add(R *res, const R *p, int pn, const R *q, int qn) {
	assert(pn >= qn);
	if(qn >= 4) {
		__m128i mod = _mm_set1_epi32(R::Mod);
		for(int i = 0; i + 3 < qn; i += 4) {
			__m128i x = _mm_loadu_si128(reinterpret_cast<const __m128i*>(p + i));
			__m128i y = _mm_loadu_si128(reinterpret_cast<const __m128i*>(q + i));
			__m128i result = add_mod_128i(x, y, mod);
			_mm_storeu_si128(reinterpret_cast<__m128i*>(res + i), result);
		}
	}
	for(int i = qn - (qn & 3); i < qn; ++ i)
		res[i] = p[i] + q[i];
	_copy(res + qn, p + qn, pn - qn);
}

void Polynomial::_subtract(R *res, const R *p, int pn, const R *q, int qn) {
	assert(pn >= qn);
	if(qn >= 4) {
		__m128i zero = _mm_setzero_si128();
		__m128i mod = _mm_set1_epi32(R::Mod);
		for(int i = 0; i + 3 < qn; i += 4) {
			__m128i x = _mm_loadu_si128(reinterpret_cast<const __m128i*>(p + i));
			__m128i y = _mm_loadu_si128(reinterpret_cast<const __m128i*>(q + i));
			__m128i result = sub_mod_128i(x, y, zero, mod);
			_mm_storeu_si128(reinterpret_cast<__m128i*>(res + i), result);
		}
	}
	for(int i = qn - (qn & 3); i < qn; ++ i)
		res[i] = p[i] - q[i];
	_copy(res + qn, p + qn, pn - qn);
}

#else

void Polynomial::_negate(R *res, const R *p, int n) {
	if(n >= 4) {
		const unsigned mods[4] = { R::Mod, R::Mod, R::Mod, R::Mod };

		__asm(
		"	vpxor %%xmm2, %%xmm2, %%xmm2;"
		"	vmovdqa %0, %%xmm3;"
		::"m"(mods));
		for(int i = 0; i + 3 < n; i += 4) {
			__asm(
			"	vpsubd	%1, %%xmm2, %%xmm0;"
			"	vpcmpgtd	%%xmm0, %%xmm2, %%xmm1;"
			"	vpand	%%xmm1, %%xmm3, %%xmm1;"
			"	vpaddd	%%xmm1, %%xmm0, %%xmm0;"
			"	vmovups	%%xmm0, %0;"
			:"=m"(res[i]):"m"(p[i]));
		}
	}
	for(int i = n - (n & 3); i < n; ++ i)
		res[i] = -p[i];
}


inline __m128i add_mod_128i(const __m128i &x, const __m128i &y, const __m128i &mod) {
	__m128i sum = _mm_add_epi32(x, y);
	__m128i lt_mod_mask = _mm_cmpgt_epi32(mod, sum);
	__m128i mod_sub = _mm_andnot_si128(lt_mod_mask, mod);
	__m128i result = _mm_sub_epi32(sum, mod_sub);
	return result;
}

inline __m128i sub_mod_128i(const __m128i &x, const __m128i &y, const __m128i &zero, const __m128i &mod) {
	__m128i diff = _mm_sub_epi32(x, y);
	__m128i lt_0_mask = _mm_cmpgt_epi32(zero, diff);
	__m128i mod_add = _mm_and_si128(lt_0_mask, mod);
	__m128i result = _mm_add_epi32(diff, mod_add);
	return result;
}
void Polynomial::_add(R *res, const R *p, int pn, const R *q, int qn) {
	assert(pn >= qn);
	if(qn >= 4) {
		const unsigned mods[4] = { R::Mod, R::Mod, R::Mod, R::Mod };
		__asm(
		"	vmovdqa %0, %%xmm2;"
		::"m"(mods));
		for(int i = 0; i + 3 < qn; i += 4) {
			__asm(
			"	vmovdqu	%1, %%xmm0;"
			"	vpaddd	%2, %%xmm0, %%xmm0;"
			"	vpcmpgtd	%%xmm0, %%xmm2, %%xmm1;"
			"	vpandn	%%xmm2, %%xmm1, %%xmm1;"
			"	vpsubd	%%xmm1, %%xmm0, %%xmm0;"
			"	vmovups	%%xmm0, %0;"
			:"=m"(res[i]):"m"(p[i]),"m"(q[i]));
		}
	}
	for(int i = qn - (qn & 3); i < qn; ++ i)
		res[i] = p[i] + q[i];
	_copy(res + qn, p + qn, pn - qn);
}

void Polynomial::_subtract(R *res, const R *p, int pn, const R *q, int qn) {
	assert(pn >= qn);
	if(qn >= 4) {
		const unsigned mods[4] = { R::Mod, R::Mod, R::Mod, R::Mod };
		__asm(
		"	vpxor %%xmm2, %%xmm2, %%xmm2;"
		"	vmovdqa %0, %%xmm3;"
		::"m"(mods));
		for(int i = 0; i + 3 < qn; i += 4) {
			__asm(
			"	vmovdqu	%1, %%xmm0;"
			"	vpsubd	%2, %%xmm0, %%xmm0;"
			"	vpcmpgtd	%%xmm0, %%xmm2, %%xmm1;"
			"	vpand	%%xmm1, %%xmm3, %%xmm1;"
			"	vpaddd	%%xmm1, %%xmm0, %%xmm0;"
			"	vmovups	%%xmm0, %0;"
			:"=m"(res[i]):"m"(p[i]),"m"(q[i]));
		}
	}
	for(int i = qn - (qn & 3); i < qn; ++ i)
		res[i] = p[i] - q[i];
	_copy(res + qn, p + qn, pn - qn);
}

#endif

void Polynomial::_add(R *p, const R *q, int n) {
	_add(p, p, n, q, n);
}

void Polynomial::_subtract(R *p, const R *q, int n) {
	_subtract(p, p, n, q, n);
}

void Polynomial::_multiply_1(R *res, const R *p, int n, R c0) {
//	assert(res != p);
	for(int i = 0; i < n; ++ i)
		res[i] = p[i] * c0;
}

void Polynomial::_multiply_1(R *p, int n, R c0) {
	_multiply_1(p, p, n, c0);
}

void Polynomial::_add_multiply_1(R *p, const R *q, int n, R c0) {
	for(int i = 0; i < n; ++ i)
		p[i] += q[i] * c0;
}

void Polynomial::_multiply_1_int(R *res, const R *p, int n, int z) {
	assert(z > 1);
	switch(z) {
	case 2: _multiply_power_of_two(res, p, n, 1); break;
	case 4: _multiply_power_of_two(res, p, n, 2); break;
	case 8: _multiply_power_of_two(res, p, n, 3); break;
	default:
		_multiply_1(res, p, n, R(z));
		break;
	}
}

void Polynomial::_divide_1_int(R *res, const R *p, int n, int z) {
	assert(z > 1);
	switch(z) {
	case 2: _divide_power_of_two(res, p, n, 1); break;
	case 4: _divide_power_of_two(res, p, n, 2); break;
	case 8: _divide_power_of_two(res, p, n, 3); break;
	default:
		_multiply_1(res, p, n, R(z).inverse());
		break;
	}
}

void Polynomial::_multiply_power_of_two(R *res, const R *p, int n, int k) {
	assert(0 < k && k < 31);
	R mul = R(1 << k);
	_multiply_1(res, p, n, mul);
}

void Polynomial::_divide_power_of_two(R *res, const R *p, int n, int k) {
	assert(0 < k && k < 31);
	static const R Inv2 = R(2).inverse();
	R inv = k == 1 ? Inv2 : R(1 << k).inverse();
	_multiply_1(res, p, n, inv);
}


void Polynomial::_multiply_select_method(R *res, const R *p, int pn, const R *q, int qn) {
	if(pn < qn) std::swap(p, q), std::swap(pn, qn);
	assert(res != p && res != q && pn >= qn && qn > 0);
	int rn = pn + qn - 1;
	if(qn == 1) {
		_multiply_1(res, p, pn, q[0]);
	}else if(qn < TOOM22_THRESHOLD) {
		_schoolbook_multiplication(res, p, pn, q, qn);
	}else if(rn < SCHONHAGE_STRASSEN_THRESHOLD) {
		int t = pn < TOOM33_THRESHOLD ? 2 : pn < TOOM44_THRESHOLD ? 3 : 4;
		int m = (pn + t - 1) / t * (t-1) + 1;
		if(qn >= m) {
			_toom44(res, p, pn, q, qn);
		}else {
			WorkSpace ws(m + (m + pn - 1));
			_copy(ws.get(), q, qn);
			_fill_zero(ws.get() + qn, m - qn);
			_toom44(ws.get() + m, p, pn, ws.get(), m);
			_copy(res, ws.get() + m, pn + qn - 1);
		}
	}else {
		_schonhage_strassen(res, p, pn, q, qn);
	}
}

void Polynomial::_square_select_method(R *res, const R *p, int pn) {
	_multiply_select_method(res, p, pn, p, pn);
}

int Polynomial::_schonhage_strassen_determine_first_k(int n) {
	static const int thresholds[][2] = {
		{ 1, 1 }, { 5, 2 }, { 25, 3 }, { 33, 2 },
		{ 49, 3 }, { 65, 2 }, { 73, 3 }, { 97, 4 },
		{ 129, 3 }, { 193, 4 }, { 257, 3 }, { 289, 4 },
		{ 385, 5 }, { 513, 4 }, { 769, 5 }, { 1025, 4 },
		{ 1153, 5 }, { 2001, 6 }, { 2049, 5 }, { 3073, 6 },
		{ 4097, 5 }, { 4609, 6 }, { 6145, 7 }, { 8193, 6 },
		{ 12289, 7 }, { 24577, 8 }, { 32769, 7 }, { 40961, 8 },
		{ 65537, 9 }, { 131073, 8 }, { 229377, 9 }, { 393217, 10 },
		{ 524289, 9 }, { 655361, 10 }, { 1048577, 11 }, { 2097153, 10 },
		{ 3670017, 11 },
	};
	static const int NumThreshold = sizeof(thresholds) / sizeof(*thresholds);

	if(n < SCHONHAGE_STRASSEN_THRESHOLD) {
		return 0;
	}else {
		int i = 0;
		while(i < NumThreshold - 1 && thresholds[i+1][0] <= n)
			++ i;
		return thresholds[i][1];
	}
}

int Polynomial::_schonhage_strassen_determine_recurse_k(int n, int k) {
	return _schonhage_strassen_determine_first_k(n);
}

#if 0
void Polynomial::_schoolbook_multiplication(R *res, const R *p, int pn, const R *q, int qn) {
	if(qn == 1) {
		_multiply_1(res, p, pn, q[0]);
		return;
	}
	assert(res != p && res != q && pn >= qn && qn > 0);

	_fill_zero(res, pn + qn - 1);
	for(int i = 0; i < pn; ++ i)
		for(int j = 0; j < qn; ++ j)
			res[i + j] += p[i] * q[j];
}
#else

inline unsigned long long emulu(unsigned x, unsigned y) {
#if defined(_MSC_VER) && !defined(MP_64)
	return __emulu(x, y);
#else
	return (unsigned long long)x * y;
#endif
}

inline unsigned getpart_epi32(__m128i x, const int i) {
	unsigned t[4];
	memcpy(t, &x, 32);
	return t[i];
}

inline unsigned long long mint_max_prod() {
	return (unsigned long long)(mint::Mod-1) * (mint::Mod-1);
}
inline unsigned long long mint_max_value() {
	return ~0ULL - mint_max_prod();
}
//mint_sub_value > ~0ULL / 2
inline unsigned long long mint_sub_value() {
	return ((~0ULL / 2) / mint::Mod + 1) * mint::Mod;
}
inline unsigned long long montgomery_max_value() {
	return ((unsigned long long)mint::Mod << mint::W) - 1;
}
inline unsigned long long montgomery_sub_value() {
	return (((montgomery_max_value() + mint_max_prod()) / 2) / mint::Mod + 1) * mint::Mod;
}

inline void schoolbook_sum_and_reduce(mint *res, const mint *p, int pn, const mint *q, int qn) {
	const int first_count = (int)min(999ULL, mint_max_value() / mint_max_prod() + 1);
	const int second_count = (int)min(999ULL, ((mint_max_value() - (mint_sub_value() - 1)) / mint_max_prod() + 1));

	for(int i = 0; i < pn + qn - 1; ++ i) {
		int jlower = i - qn + 1 >= 0 ? i - (qn - 1) : 0;
		int jupper = pn <= i ? pn - 1 : i;

#ifdef MP_64
		unsigned long long sum = 0;
		int j = jlower;
		int firstu = j + first_count;
		if(firstu >= jupper)
			goto rem_j;
		for( ; j < firstu; ++ j)
			sum += (unsigned long long)p[j].residue * q[i-j].residue;
		if(sum >= mint_sub_value())
			sum -= mint_sub_value();

		while(1){
			int u = j + second_count;
			if(u >= jupper) break;
			for( ; j < u; ++ j)
				sum += (unsigned long long)p[j].residue * q[i-j].residue;
			if(sum >= mint_sub_value())
				sum -= mint_sub_value();
		}

	rem_j:
		for(; j <= jupper; ++ j)
			sum += (unsigned long long)p[j].residue * q[i-j].residue;

		res[i].residue = mint::reduce((unsigned)(sum % mint::Mod));
#else
		unsigned long long sum = 0;

		for(int j = jlower; j <= jupper; ++ j) {
			sum += (unsigned long long)p[j].residue * q[i-j].residue;
			if(sum >= montgomery_sub_value())
				sum -= montgomery_sub_value();
		}
		res[i].residue = mint::reduce(sum);
#endif
	}

}

void Polynomial::_schoolbook_multiplication(R *res, const R *p, int pn, const R *q, int qn) {
	assert(res != p && res != q && pn >= qn && qn > 1);
	schoolbook_sum_and_reduce(res, p, pn, q, qn);
}

#endif

void Polynomial::_toom22(R *res, const R *p, int pn, const R *q, int qn, WorkSpaceStack &wss) {
	assert(res != p && res != q && pn >= qn);
	assert(qn >= TOOM22_THRESHOLD);

	int lo = (pn + 1) / 2, hp = pn - lo, hq = qn - lo;
	assert(0 < lo && 0 < hq && hp <= lo && hq <= lo);
	//r0 = pq(0), r1 = pq(1), rinf = pq(inf)
	//w0 = r0, w1 = -r0 + r1 - rinf, w2 = rinf
	WorkSpace ws(wss, lo * 4 - 1);
	R *t0 = ws.get(), *t1 = t0 + lo;
	R *r0 = res, *r1 = t1 + lo, *rinf = res + lo * 2;
	_add(t0, p, lo, p + lo, hp);
	_add(t1, q, lo, q + lo, hq);
	if(lo < TOOM22_THRESHOLD) {
		_schoolbook_multiplication(r1, t0, lo, t1, lo);
		_schoolbook_multiplication(r0, p, lo, q, lo);
	}else {
		_toom22(r1, t0, lo, t1, lo, wss);
		_toom22(r0, p, lo, q, lo, wss);
	}
	res[lo + lo - 1] = ZeroR();
	if(hq < TOOM22_THRESHOLD) {
		_schoolbook_multiplication(rinf, p + lo, hp, q + lo, hq);
	}else {
		_toom22(rinf, p + lo, hp, q + lo, hq, wss);
	}
	_subtract(r1, r0, lo + lo - 1);
	_subtract(r1, rinf, hp + hq - 1);
	_add(res + lo, r1, lo + lo - 1);
}

int Polynomial::_toom22_estimate_workspace(int pn, int qn) {
	assert(pn >= qn);
	if(qn < TOOM22_THRESHOLD)
		return 0;
	int lo = (pn + 1) / 2;
	return lo * 4 - 1 + _toom22_estimate_workspace(lo, lo);
}

void Polynomial::_toom33(R *res, const R *p, int pn, const R *q, int qn) {
	assert(res != p && res != q && pn >= qn);

	if(pn < TOOM33_THRESHOLD) {
		if(pn < TOOM22_THRESHOLD) {
			_schoolbook_multiplication(res, p, pn, q, qn);
		}else {
			WorkSpaceStack wss(_toom22_estimate_workspace(pn, qn));
			_toom22(res, p, pn, q, qn, wss);
		}
		return;
	}

	//r0   = pq(0)   =  p0                *  q0
	//r1   = pq(1)   = (p0 +   p1 +   p2) * (q0 +   q1 +   q2)
	//rm1  = pq(-1)  = (p0 -   p1 +   p2) * (q0 -   q1 +   q2)
	//r2   = pq(2)   = (p0 + 2*p1 + 4*p2) * (q0 + 2*q1 + 4*q2)
	//rinf = pq(inf) =                p2  *                q2

	//w0 = r0
	//w1 = r1 + 2*rinf - ((2*rm1 + r2)/3 + r0)/2
	//w2 = (r1 + rm1)/2 - r0 - rinf
	//w3 = (r0 - r1 + (r2 - rm1)/3)/2 - 2*rinf
	//w4 = rinf

	int lo = (pn + 2) / 3, hp = pn - lo * 2, hq = qn - lo * 2;
	assert(0 < lo && 0 < hq && hp <= lo && hq <= lo);

	Polynomial p0(p, p + lo), p1(p + lo, p + lo * 2), p2(p + lo * 2, p + pn);
	Polynomial q0(q, q + lo), q1(q + lo, q + lo * 2), q2(q + lo * 2, q + qn);

	Polynomial r0 = p0 * q0;
	Polynomial r1 = (p0 + p1 + p2) * (q0 + q1 + q2);
	Polynomial rm1 = (p0 - p1 + p2) * (q0 - q1 + q2);
	Polynomial r2 = (p0 + p1*2 + p2*4) * (q0 + q1*2 + q2*4);
	Polynomial rinf = p2 * q2;

	static const R R_1_2 = R(2).inverse(), R_1_3 = R(3).inverse(), R_1_6 = R_1_2 * R_1_3;

	Polynomial r0_1_2 = r0 * R_1_2;
	Polynomial r1_1_2 = r1 * R_1_2;
	Polynomial r2_1_6 = r2 * R_1_6;
	Polynomial rinf_2 = rinf * 2;

	Polynomial w0 = r0;
	Polynomial w1 = r1 - r0_1_2 - rm1 * R_1_3 - r2_1_6 + rinf_2;
	Polynomial w2 = r1_1_2 - r0 + rm1 * R_1_2 - rinf;
	Polynomial w3 = r0_1_2 - r1_1_2 - rm1 * R_1_6 + r2_1_6 - rinf_2;
	Polynomial w4 = rinf;
	w3.resize(lo + hp - 1);

	Polynomial w = w0 + (w1 << lo) + (w2 << (lo * 2)) + (w3 << (lo * 3)) + (w4 << (lo * 4));

	assert(w.size() == pn + qn - 1);
	_copy(res, w.data(), pn + qn - 1);
}

void Polynomial::_toom44(R *res, const R *p, int pn, const R *q, int qn) {
	assert(res != p && res != q && pn >= qn);

	if(pn < TOOM44_THRESHOLD) {
		_toom33(res, p, pn, q, qn);
		return;
	}

	int lo = (pn + 3) / 4, hp = pn - lo * 3, hq = qn - lo * 3;
	assert(0 < lo && 0 < hq && hp <= lo && hq <= lo);

	Polynomial p0(p, p + lo), p1(p + lo, p + lo * 2), p2(p + lo * 2, p + lo * 3), p3(p + lo * 3, p + pn);
	Polynomial q0(q, q + lo), q1(q + lo, q + lo * 2), q2(q + lo * 2, q + lo * 3), q3(q + lo * 3, q + qn);

	static const R R_1_2 = R(2).inverse(), R_1_3 = R(3).inverse(), R_1_9 = R_1_3 * R_1_3;
	static const R R_1_4 = R(4).inverse(), R_1_6 = R(6).inverse(), R_1_18 = R_1_9 * R_1_2, R_1_45 = R(45).inverse(), R_1_90 = R(90).inverse();
	static const R R_3_2 = R(3) / 2, R_5_2 = R(5) / 2, R_5_4 = R_5_2 * R_1_2;
	static const R R_7_18 = R(7) / 18, R_32_9 = R(32) / 9, R_2_3 = R(2) / 3, R_2_9 = R(2) / 9, R_2_45 = R(2) / 45;
	static const R R_16_15 = R(16) / 15, R_8_3 = R(8) / 3, R_16_9 = R(16) / 9;

	Polynomial r0 = p0 * q0;
	Polynomial rh = (p0 + (p1 + (p2 + p3/2)/2)/2) * (q0 + (q1 + (q2 + q3/2)/2)/2);
	Polynomial rmh = (p0 + (-p1 + (p2 - p3/2)/2)/2) * (q0 + (-q1 + (q2 - q3/2)/2)/2);
	Polynomial r1 = (p0 + p1 + p2 + p3) * (q0 + q1 + q2 + q3);
	Polynomial rm1 = (p0 - p1 + p2 - p3) * (q0 - q1 + q2 - q3);
	Polynomial r2 = (p0 + (p1 + (p2 + p3*2)*2)*2) * (q0 + (q1 + (q2 + q3*2)*2)*2);
	Polynomial rinf = p3 * q3;

	Polynomial r1_2_3 = r1 * R_2_3;
	Polynomial rmh_16_15 = rmh * R_16_15;
	Polynomial rmh_8_3 = rmh * R_8_3;
	Polynomial rh_16_9 = rh * R_16_9;
	Polynomial rh_8_3 = rh * R_8_3;

	Polynomial w0 = r0;
	Polynomial w1 = -r0*R_1_2 - r1*R_1_3 - rinf*R_1_2 + rm1*R_1_9 + r2*R_1_90 + rh_16_9 - rmh_16_15;
	Polynomial w2 = -r0*5 - r1*R_1_6 + rinf*R_1_4 - rm1*R_1_6 + rh_8_3 + rmh_8_3;
	Polynomial w3 = r0*R_5_2 + r1*R_3_2 + rinf*R_5_2 - rm1*R_7_18 - r2*R_1_18 - rh*R_32_9;
	Polynomial w4 = r0*4 + r1_2_3 - rinf*R_5_4 + rm1*R_2_3 - rh_8_3 - rmh_8_3;
	Polynomial w5 = -r0*2 - r1_2_3 - rinf*2 - rm1*R_2_9 + r2*R_2_45 + rh_16_9 + rmh_16_15;
	Polynomial w6 = rinf;
	w5.resize(lo + hp - 1);

	Polynomial w = w0 + (w1 << lo) + (w2 << (lo * 2)) + (w3 << (lo * 3)) + (w4 << (lo * 4)) + (w5 << (lo * 5)) + (w6 << (lo * 6));

	assert(w.size() == pn + qn - 1);
	_copy(res, w.data(), pn + qn - 1);
}

void Polynomial::_schonhage_strassen(R *res, const R *p, int pn, const R *q, int qn) {
	int rn = pn + qn - 1;
	int k = _schonhage_strassen_determine_first_k(rn);
	assert(k > 0);
	int K = 1 << k;
	int N = (rn + K - 1) >> k << k;

	WorkSpace ws(N * 3);
	R *Z = ws.get(), *X = Z + N, *Y = X + N;

	_copy(X, p, pn);
	_fill_zero(X + pn, N - pn);
	_copy(Y, q, qn);
	_fill_zero(Y + qn, N - qn);

	_schonhage_strassen_recurse(Z, X, Y, N, k);

	_copy(res, Z, rn);
}

void Polynomial::_schonhage_strassen_recurse(R *res, const R *X, const R *Y, int N, int k) {
	assert(1 < k && k <= 15);
	assert(0 < N && (N & ((1 << k)-1)) == 0);
	int K = 1 << k;
	int m = N >> k;
	int n = (m + m - 1 + K - 1) >> k << k, n2 = n * 2;
	int nK = n << k;
	assert(n < N);
//	cerr << "s.s. recurse N = " << N << ", k = " << k << ", n = " << n << endl;

	WorkSpace ws(nK * 2 + n * 2, true);
	R *x = ws.get(), *y = x + nK, *tmp = y + nK;

	_schonhage_strassen_decompose(x, X, k, n, m, tmp);
	_schonhage_strassen_decompose(y, Y, k, n, m, tmp);
	
	unsigned short *omega_table = new unsigned short[K - 1];
	_fft_make_omega_table(n, k, omega_table);

	_fft_rec(x, n, n, k, tmp, omega_table);
	_fft_rec(y, n, n, k, tmp, omega_table);

	int nextk = _schonhage_strassen_determine_recurse_k(n, k);
	assert(nextk <= k);

	if(nextk == 0) {
		for(int i = 0; i < nK; i += n) {
			_toom44(tmp, x + i, n, y + i, n);
			_copy(x + i, tmp, n);
			_subtract(x + i, tmp + n, n - 1);
		}
	}else {
		for(int i = 0; i < nK; i += n) {
			_schonhage_strassen_recurse(x + i, x + i, y + i, n, nextk);
		}
	}

	_fft_inverse_rec(x, n, k, tmp);
	_divide_power_of_two(x, x, nK, k);

	delete[] omega_table;

	int inv_omega_k = n2 - (n >> k);

	_fill_zero(res, N);
	for(int i = 0, j = 0, shift = 0; i < nK; i += n, j += m) {
		if(j >= 2 * N) j -= 2 * N;
		_negacyclic_shift(tmp, x + i, n, shift);
		_negacyclic_shift_add(res, N, tmp, n, j);
		shift += inv_omega_k;
		if(shift >= n2) shift -= n2;
	}
}

void Polynomial::_schonhage_strassen_decompose(R *x, const R *X, int k, int n, int m, R *tmp) {
	int N = m << k, n2 = n * 2, omega = n >> k;
	for(int i = 0, j = 0, shift = 0; j < N; i += n, j += m) {
		_copy(tmp, X + j, m);
		_fill_zero(tmp + m, n - m);
		_negacyclic_shift(x + i, tmp, n, shift);
		shift += omega;
		if(shift >= n2) shift -= n2;
	}
}

void Polynomial::_fft_make_omega_table(int n, int k, unsigned short *omega_table) {
	int n2 = n * 2;
	unsigned short *row = omega_table;
	for(int i = k; i > 0; -- i) {
		int K = 1 << i, hK = 1 << (i-1), hhK = hK >> 1;
		int omega = n2 >> i;
		int r = 0;
		for(int j = 0; j < hK; ++ j) {
			//r = bitrev(j * 2, i) = bitrev(j, i-1)
			row[j] = r * omega % n2;
			for(int h = hhK; h > (r ^= h); h >>= 1) ;
		}
		row += hK;
	}
}

void Polynomial::_fft_rec(R *x, int n, int inc, int k, R *tmp, const unsigned short *omega_table) {
	if(k == 0) return;
	int hK = 1 << (k-1);
	R *x0 = x, *x1 = x + inc;
	_fft_rec(x0, n, inc * 2, k - 1, tmp, omega_table + hK);
	_fft_rec(x1, n, inc * 2, k - 1, tmp, omega_table + hK);

	for(int i = 0; i < hK; ++ i, x0 += inc * 2, x1 += inc * 2) {
		int shift = omega_table[i];
		_negacyclic_shift(tmp, x1, n, shift);
		_subtract(x1, x0, n, tmp, n);
		_add(x0, tmp, n);
	}
}

void Polynomial::_fft_inverse_rec(R *x, int n, int k, R *tmp) {
	if(k == 0) return;
	int hK = 1 << (k-1);
	R *x0 = x, *x1 = x + (n << (k-1));
	_fft_inverse_rec(x0, n, k - 1, tmp);
	_fft_inverse_rec(x1, n, k - 1, tmp);

	int n2 = n * 2, omega = n2 - (n2 >> k);
	for(int i = 0, shift = 0; i < hK; ++ i, x0 += n, x1 += n) {
		_negacyclic_shift(tmp, x1, n, shift);
		_subtract(x1, x0, n, tmp, n);
		_add(x0, tmp, n);
		shift += omega;
		if(shift >= n2) shift -= n2;
	}
}

void Polynomial::_negacyclic_shift(R *res, const R *x, int n, int shift) {
	assert(0 <= shift && shift <= n * 2);
	if(shift < n) {
		_negate(res, x + (n - shift), shift);
		_copy(res + shift, x, n - shift);
	}else {
		int s = shift - n;
		_copy(res, x + (n - s), s);
		_negate(res + s, x, n - s);
	}
}

void Polynomial::_negacyclic_shift_add(R *p, int pn, const R *q, int qn, int shift) {
	assert(p != q && pn >= qn && 0 <= shift && shift < pn * 2);
	if(shift < pn) {
		if(shift + qn <= pn) {
			_add(p + shift, q, qn);
		}else {
			int t = pn - shift;
			_subtract(p, q + t, qn - t);
			_add(p + shift, q, t);
		}
	}else {
		int s = shift - pn;
		if(s + qn <= pn) {
			_subtract(p + s, q, qn);
		}else {
			int t = pn - s;
			_add(p, q + t, qn - t);
			_subtract(p + s, q, t);
		}
	}
}


void Polynomial::_reverse(R *res, const R *p, int pn) {
	if(res == p) {
		std::reverse(res, res + pn);
	}else {
		for(int i = 0; i < pn; ++ i)
			res[pn - 1 - i] = p[i];
	}
}

void Polynomial::_inverse_power_series(R *res, int resn, const R *p, int pn) {
	if(resn == 0) return;
	assert(res != p);
	assert(p[0] == OneR());
	R two = R(2);
	WorkSpace ws(resn * 4);
	R *tmp1 = ws.get(), *tmp2 = tmp1 + resn * 2;
	_fill_zero(res, resn);
	res[0] = p[0];
	int curn = 1;
	while(curn < resn) {
		int nextn = std::min(resn, curn * 2);
		_square_select_method(tmp1, res, curn);
		_multiply_select_method(tmp2, tmp1, std::min(nextn, curn * 2 - 1), p, std::min(nextn, pn));
		_multiply_power_of_two(res, res, curn, 1);
		_subtract(res, tmp2, nextn);
		curn = nextn;
	}
}

void Polynomial::_precompute_inverse(R *res, int resn, const R *p, int pn) {
	WorkSpace ws(pn);
	R *tmp = ws.get();
	_reverse(tmp, p, pn);
	_inverse_power_series(res, resn, tmp, pn);
}

void Polynomial::_divide_precomputed_inverse(R *res, int resn, const R *revp, int pn, const R *inv) {
	WorkSpace ws(pn + resn);
	R *tmp = ws.get();
	_multiply_select_method(tmp, revp, pn, inv, resn);
	_reverse(res, tmp, resn);
}

void Polynomial::_divide_remainder_precomputed_inverse(R *quot, R *rem, const R *p, int pn, const R *q, int qn, const R *inv) {
	if(pn < qn) {
		_copy(rem, p, pn);
		_fill_zero(rem + pn, qn - 1 - pn);
		return;
	}
	assert(qn > 0);
	assert(q[qn - 1] == OneR());
	if(qn == 1) return;
	int quotn = pn - qn + 1;
	int rn = qn - 1, tn = std::min(quotn, rn), un = tn + rn;
	WorkSpace ws(pn + un + (quot != NULL ? 0 : quotn));
	R *revp = ws.get(), *quotmul = revp + pn;
	if(quot == NULL) quot = quotmul + un;
	_reverse(revp, p, pn);
	_divide_precomputed_inverse(quot, quotn, revp, pn, inv);
	_multiply_select_method(quotmul, q, rn, quot, tn);
	_subtract(rem, p, rn, quotmul, rn);
}


}

//ループの順序を変えただけでとてつもなく速度が変わった…キャッシュか…
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;
	int z = a[0];
	rep(k, 6) {
		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.residue != 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_polynomial;
	int n = (int)c.size();
	Polynomial poly;
	poly.resize(n+1);
	rep(i, n)
		poly.set(i, -c[i]);
	poly.set(n, 1);

	Polynomial resp = Polynomial::X().powerMod(N, poly);

	vector<mint> res(n);
	rep(i, n)
		res[i] = resp.get(i);
	return res;
}

int main() {
	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, cntC;
		cntP = doDP(primes, P);
		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;
		coeffs = computeLinearRecurrentSequenceCoefficients(c, X + (N - 1));
		mint ans = accumulate(all(coeffs), mint());
		printf("%d\n", ans.get());
	}
	return 0;
}
0