#919158 (C++23) No.754 畳み込みの和

提出ソース
結果

問題	No.754 畳み込みの和
ユーザー	popofy
提出日時	2023-10-09 00:38:34
言語	C++23 (gcc 13.3.0 + boost 1.87.0)
結果	AC
実行時間	223 ms / 5,000 ms
コード長	20,425 bytes
コンパイル時間	2,845 ms
コンパイル使用メモリ	218,740 KB
実行使用メモリ	18,524 KB
最終ジャッジ日時	2024-07-26 18:18:12
合計ジャッジ時間	4,790 ms
ジャッジサーバーID （参考情報）	judge4 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン	結果
other	AC * 3
権限があれば一括ダウンロードができます
ソースコード

raw source code
// #define Q__OPTIMIZE
// #define Q__INCLUDE_ATCODER_LIB
// #define Q__INTERACTIVE
#if !__INCLUDE_LEVEL__
#include __FILE__



struct Solver {
  void solve() {
    INT(n);
    VEC(ll,a,n+1); VEC(ll,b,n+1);
    Polynomial<ll> P,Q;
    P.coef = a;
    Q.coef = b;
    auto C = P * Q;
    DUMP(C.coef);
    ll ans = 0;
    REP(i,n+1) ans = safe_mod(ans+C[i],MOD);
    print(ans);
  }


  void naive() {

  }
} solver;



signed main(void){
  NO_SYNC_STD;
  V<string> options;
#ifdef Q__OPTIMIZE
  options.push_back("OPTIMIZE");
#endif
#ifdef Q__INTERACTIVE
  options.push_back("INTERACTIVE");
#endif
#ifdef Q__INCLUDE_ATCODER_LIB
  options.push_back("INCLUDE_ATCODER_LIB");
#endif
  DUMP(options);
#ifndef Q__NAIVE
  solver.solve();
#else
  DUMP("naive");
  solver.naive();
#endif
  return 0;
}

#else
#define _GLIBCXX_DEQUE_BUF_SIZE 64
#ifdef Q__OPTIMIZE
#pragma GCC target("avx")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")
#endif
#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <deque>
#include <forward_list>
#include <fstream>
#include <functional>
#include <iomanip>
#include <ios>
#include <iostream>
#include <limits>
#include <list>
#include <map>
#include <memory>
#include <numeric>
#include <optional>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <thread>
#include <tuple>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#ifdef Q__INCLUDE_ATCODER_LIB
#include <atcoder/all>
using namespace atcoder;
//using mint = modint1000000007;
using mint = modint998244353;
std::istream &operator>>(std::istream& is, mint& a) { long long tmp; is >> tmp; a = tmp; return is; }
std::ostream &operator<<(std::ostream& os, const mint& a) {os << a.val(); return os;}
#endif
using namespace std;
#define MOD 1000000007
#define OVERLOAD4(a, b, c, d, e, ...) e
#define REP1(a)          for(decltype(a) i = 0, i##_len = (a); i < i##_len; ++i)
#define REP2(i, a)       for(decltype(a) i = 0, i##_len = (a); i < i##_len; ++i)
#define REP3(i, a, b)    for(decltype(b) i = (a), i##_len = (b); i < i##_len; ++i)
#define REP4(i, a, b, c) for(decltype(b) i = (a), i##_len = (b); i < i##_len; i += (c))
#define REP(...) OVERLOAD4(__VA_ARGS__, REP4, REP3, REP2, REP1)(__VA_ARGS__)
#define RREP1(a)          for(decltype(a) i = (a); i--;)
#define RREP2(i, a)       for(decltype(a) i = (a); i--;)
#define RREP3(i, a, b)    for(decltype(a) i = (b), i##_len = (a); i-- > i##_len;)
#define RREP4(i, a, b, c) for(decltype(a) i = (a)+((b)-(a)-1)/(c)*(c), i##_len = (a); i >= i##_len; i -= c)
#define RREP(...) OVERLOAD4(__VA_ARGS__, RREP4, RREP3, RREP2, RREP1)(__VA_ARGS__)
#define MREP(v,...) for(auto v:make_enum_vec({__VA_ARGS__}))
#define QREP(q, l, r, n) for (ll q = 1, l = n / (q + 1) + 1, r = n / q + 1; q <= n; q = (q == n ? n + 1 : n / (n / (q + 1))), l = n / (q + 1) + 1, r = n / q + 1)
#define COMB_REP(i,n,k) for (ll t, i = POW2(k) - 1; i < POW2(n); t=i|(i-1), i = (t+1)|(((~t & - ~t)-1) >> (__builtin_ctzll(i)+1)))
#define SUBSET_ENUM_REP(i,s) for (ll i = (1LL << 60) - 1; i >= 0, i &= s; --i)
#define SUBSET_INCLUDE_REP(i,n,s) for (int i = s; i < POW2(n); i=(++i)|s)
#define POPONLY_REP(i,s) for (ll i=s&-s; i; i=s&(~s+(i << 1)))
#define ALL(x)  (x).begin(), (x).end()
#define RALL(x) (x).rbegin(), (x).rend()
#define SZ(x)   ((int)(x).size())
#define POW2(n)      (1LL << ((int)(n)))
#define GET1BIT(x,n) (((x) >> (int)(n)) & 1)
#define INF ((1 << 30) - 1)
#define INFL (1LL << 60)
#define PRECISION std::setprecision(16)
#define SLEEP(n) std::this_thread::sleep_for(std::chrono::seconds(n))
#define INT(...) int __VA_ARGS__;    input(__VA_ARGS__)
#define LL(...)  ll __VA_ARGS__;     input(__VA_ARGS__)
#define STR(...) string __VA_ARGS__; input(__VA_ARGS__)
#define LD(...)  ld __VA_ARGS__;     input(__VA_ARGS__)
#define VEC(type, name, size) vector<type> name(size); input(name)
#ifdef Q__INTERACTIVE
#define NO_SYNC_STD
#define ENDL std::endl
#else
#define NO_SYNC_STD std::cin.tie(nullptr);ios::sync_with_stdio(false)
#define ENDL "\n"
#endif
#ifdef Q__LOCAL
#include <dump.hpp>
#define DUMP(...) DUMPOUT << "  " << string(#__VA_ARGS__) << ": " << "[" << to_string(__LINE__) << ":" << __FUNCTION__ << "]" << endl ,dump_func(__VA_ARGS__)
#define VDUMP(...) DUMPOUT << "  " << string(#__VA_ARGS__) << ": " << "[" << to_string(__LINE__) << ":" << __FUNCTION__ << "]" << endl, vdump_func(__VA_ARGS__)
#else
#define DUMP(...)
#define VDUMP(...)
#endif
using ll=long long;
using ull=unsigned long long;
using ld=long double;
template<class T> using V=vector<T>;
template<class T> using VV=vector<vector<T>>;
template<class T> using PQ=priority_queue<T,V<T>,greater<T>>;
template<class T> istream &operator>>(istream &is,V<T> &v){for(auto&& e:v)is >> e;return is;}
template<class T> istream &operator>>(istream &is,complex<T> &v){T x,y; is >> x >> y;v.real(x);v.imag(y);return is;}
template<class T,class U> istream &operator>>(istream &is,pair<T,U> &v){is >> v.first >> v.second;return is;}
template<class T,size_t n> istream &operator>>(istream &is,array<T,n> &v){for(auto&& e:v)is >> e;return is;}
template<class... A> void input(A&&... args){(cin >> ... >> args);}
template<class... A> void print_rest(){cout << ENDL;}
template<class T,class... A> void print_rest(const T& first,const A&... rest){cout << " " << first;print_rest(rest...);}
template<class T,class... A> void print(const T& first,const A&... rest){cout << fixed << PRECISION << first;print_rest(rest...);}
template<class T,class... A> void die(const T& first,const A&... rest){cout << fixed << PRECISION << first;print_rest(rest...);exit(0);}
template <typename ... Args> string fmt(const string& fmt, Args ... args ){size_t len = snprintf( nullptr, 0, fmt.c_str(), args ... );vector<char> buf(len + 1);snprintf(&buf[0], len + 1, fmt.c_str(), args ... );return string(&buf[0], &buf[0] + len);}
template<class T> inline string join(const T& v,string sep=" "){if(!SZ(v))return "";stringstream ss;for(auto&& e:v)ss << sep << e;return ss.str().substr(SZ(sep));}
V<string> split(const string &s,char sep=' ') {V<string> ret;stringstream ss(s);string buf;while(getline(ss,buf,sep))ret.push_back(buf);return ret;}
template<class T> inline string padding(const T& v,int len,char pad=' ',bool l=false){stringstream ss;ss << (l?std::left:std::right) << setw(len) << setfill(pad) << v;return ss.str();}
template<class T> V<T> make_vec(size_t n,T a){return V<T>(n,a);}
template<class... Ts> auto make_vec(size_t n,Ts... ts){return V<decltype(make_vec(ts...))>(n,make_vec(ts...));}
template<class T> inline bool chmax(T& a,T b){if(a<b){a=b;return 1;} return 0;}
template<class T> inline bool chmin(T& a,T b){if(a>b){a=b;return 1;} return 0;}
template<class T,class F> pair<T,T> binarysearch(T ng,T ok,T eps,F f,bool sign=false){while(abs(ng-ok)>eps){auto mid=ng+(ok-ng)/2;if(sign^f(mid)){ok=mid;}else{ng=mid;}}return{ng,ok};}
template<class T> constexpr T cdiv(T x,T y){return (x+y-1)/y;}
template<class T> constexpr bool between(T a,T x,T b){return(a<=x&&x<b);}
template<class T> constexpr T pos1d(T y,T x,T h,T w){assert(between(T(0),y,h));assert(between(T(0),x,w));return y*w+x;}
template<class T> constexpr pair<T,T> pos2d(T p,T h,T w){T y=p/w,x=p-y*w;assert(between(T(0),y,h));assert(between(T(0),x,w));return{y,x};}
template<class T> constexpr T sign(T n) {return (n > 0) - (n < 0);}
template<class T> inline V<T> transposed(V<T>& A){int h=SZ(A),w=SZ(A[0]);V<T> tA(w);REP(i,h)REP(j,w)tA[j].push_back(A[i][j]);return tA;}
template<class T> inline V<T> ruiseki(V<T>& a){auto ret = a; ret.push_back(T(0));exclusive_scan(ALL(ret), ret.begin(), 0);return ret;}
template<class T> inline V<T> kaisa(V<T>& a){V<T> ret(a.size());adjacent_difference(ALL(a), ret.begin());return ret;}
template<class T> inline int g_index(V<T> &s, T x) {
  if (s.empty()) return -2;
  auto it = upper_bound(ALL(s), x);
  if (it == s.end()) return -1;
  return (int)distance(s.begin(), it);
}
template<class T> inline int ge_index(V<T> &s, T x) {
  if (s.empty()) return -2;
  auto it = lower_bound(ALL(s), x);
  if (it == s.end()) return -1;
  return (int)distance(s.begin(), it);
}
template<class T> inline int l_index(V<T> &s, T x) {
  if (s.empty()) return -2;
  auto it = lower_bound(ALL(s), x);
  if (it == s.begin()) return -1;
  return (int)distance(s.begin(), prev(it));
}
template<class T> inline int le_index(V<T> &s, T x) {
  if (s.empty()) return -2;
  auto it = upper_bound(ALL(s), x);
  if (it == s.begin()) return -1;
  return (int)distance(s.begin(), prev(it));
}
template<class T> inline pair<typename set<T>::iterator,bool> g_it(set<T> &s, T x) {
  if (s.empty()) return {s.end(), false};
  auto it = s.upper_bound(x);
  if (it == s.end()) return {it, false};
  return {it, true};
}
template<class T> inline pair<typename set<T>::iterator,bool> ge_it(set<T> &s, T x) {
  if (s.empty()) return {s.end(), false};
  auto it = s.lower_bound(x);
  if (it == s.end()) return {it, false};
  return {it, true};
}
template<class T> inline pair<typename set<T>::iterator,bool> l_it(set<T> &s, T x) {
  if (s.empty()) return {s.end(), false};
  auto it = s.lower_bound(x);
  if (it == s.begin()) return {it, false};
  return {prev(it), true};
}
template<class T> inline pair<typename set<T>::iterator,bool> le_it(set<T> &s, T x) {
  if (s.empty()) return {s.end(), false};
  auto it = s.upper_bound(x);
  if (it == s.begin()) return {it, false};
  return {prev(it), true};
}
template<class T> inline V<T> it_range(set<T> &st, int l, int r) {
  auto startIt = st.lower_bound(l); auto endIt = st.upper_bound(r); V<T> ret;
  for(auto itr = startIt; itr != endIt; itr++) ret.emplace_back(*itr);
  return ret;
}
template<class T> inline pair<typename multiset<T>::iterator,bool> g_it(multiset<T> &s, T x) {
  if (s.empty()) return {s.end(), false};
  auto it = s.upper_bound(x);
  if (it == s.end()) return {it, false};
  return {it, true};
}
template<class T> inline pair<typename multiset<T>::iterator,bool> ge_it(multiset<T> &s, T x) {
  if (s.empty()) return {s.end(), false};
  auto it = s.lower_bound(x);
  if (it == s.end()) return {it, false};
  return {it, true};
}
template<class T> inline pair<typename multiset<T>::iterator,bool> l_it(multiset<T> &s, T x) {
  if (s.empty()) return {s.end(), false};
  auto it = s.lower_bound(x);
  if (it == s.begin()) return {it, false};
  return {prev(it), true};
}
template<class T> inline pair<typename multiset<T>::iterator,bool> le_it(multiset<T> &s, T x) {
  if (s.empty()) return {s.end(), false};
  auto it = s.upper_bound(x);
  if (it == s.begin()) return {it, false};
  return {prev(it), true};
}
template<class T> inline V<T> it_range(multiset<T> &st, int l, int r) {
  auto startIt = st.lower_bound(l); auto endIt = st.upper_bound(r); V<T> ret;
  for(auto itr = startIt; itr != endIt; itr++) ret.emplace_back(*itr);
  return ret;
}
constexpr ll modpow(ll x,ll n,ll m=1152921504606846976LL){ll ret=1;for(;n>0;x=x*x%m,n>>=1)if(n&1)ret=ret*x%m;return ret;}
constexpr ll safe_mod(ll x, ll m) {x%=m;if(x<0)x+=m;return x;}
constexpr ll keta(ll n, ll base = 10LL) {ll ret = 0; while(n > 0) {n /= base, ret++;} return ret;}
constexpr int pcnt(ll x) {return __builtin_popcountll(x);}
constexpr int log2f(ll x) {return 63 - __builtin_clzll(x);}
constexpr int log2c(ll x) {return (x==1LL)?0:(64-__builtin_clzll(x-1LL));}
constexpr ll nC2(ll n) {return n*(n-1)/2;}
constexpr ld deg2rad(ll degree){return (ld)degree * M_PI/180;}
mt19937 rnd_engine{random_device{}()};
inline int rand(int l, int r) {uniform_int_distribution<> ret(l, r);return ret(rnd_engine);}
inline ld lrand(ld l, ld r) {uniform_real_distribution<> ret(l, r);return ret(rnd_engine);}
inline ld nrand(ld ave, ld var) {normal_distribution<> ret(ave, var);return ret(rnd_engine);}
inline void yes(bool cond) {cout << (cond?"Yes":"No") << ENDL;}
inline bool is_palindrome(const string& s){return equal(ALL(s), s.rbegin());}
inline string make_palindrome(const string& s, bool odd = true) {string t = s.substr(0, SZ(s)-odd);reverse(ALL(t));return s + t;}
VV<int> make_enum_vec(V<int> v){
  if(v.empty()) return VV<int>(1,V<int>());
  int n=v.back(); v.pop_back();
  VV<int> ret,tmp=make_enum_vec(v);
  for(auto e:tmp)for(int i=0;i<n;++i){ret.push_back(e);ret.back().push_back(i);}
  return ret;
}
V<int> restore_path(V<int>& to, int goal, bool to1indexed = true) {
  V<int> ret;
  int x = goal;
  while(x >= 0) {
    ret.push_back(x);
    x = to[x];
  }
  reverse(ALL(ret));
  if (to1indexed) for(auto&& e: ret) e++;
  return ret;
}
const int dx4[4] = {1, 0, -1, 0};
const int dy4[4] = {0, 1, 0, -1};
const int dx6[6] = {1, 0, -1, 0, 1, -1};
const int dy6[6] = {0, 1, 0, -1, 1, -1};
const int dx8[8] = {1, 0, -1, 0, 1, -1, -1, 1};
const int dy8[8] = {0, 1, 0, -1, 1, 1, -1, -1};



template<class T> T extgcd(T a, T b, T& x, T& y) { for (T u = y = 1, v = x = 0; a;) { T q = b / a; swap(x -= q * u, u); swap(y -= q * v, v); swap(b -= q * a, a); } return b; }
template<class T> T mod_inv(T a, T m) { T x, y; extgcd(a, m, x, y); return (m + x % m) % m; }

template<int mod, int primitive_root>
class NTT {
public:
	int get_mod() const { return mod; }
	void _ntt(vector<ll>& a, int sign) {
		const int n = SZ(a);
		assert((n ^ (n&-n)) == 0); //n = 2^k

		const int g = 3; //g is primitive root of mod
		int h = (int)modpow(g, (mod - 1) / n, mod); // h^n = 1
		if (sign == -1) h = (int)mod_inv(h, mod); //h = h^-1 % mod

		//bit reverse
		int i = 0;
		for (int j = 1; j < n - 1; ++j) {
			for (int k = n >> 1; k >(i ^= k); k >>= 1);
			if (j < i) swap(a[i], a[j]);
		}

		for (int m = 1; m < n; m *= 2) {
			const int m2 = 2 * m;
			const ll base = modpow(h, n / m2, mod);
			ll w = 1;
			REP(x, m) {
				for (int s = x; s < n; s += m2) {
					ll u = a[s];
					ll d = a[s + m] * w % mod;
					a[s] = u + d;
					if (a[s] >= mod) a[s] -= mod;
					a[s + m] = u - d;
					if (a[s + m] < 0) a[s + m] += mod;
				}
				w = w * base % mod;
			}
		}

		for (auto& x : a) if (x < 0) x += mod;
	}
	void ntt(vector<ll>& input) {
		_ntt(input, 1);
	}
	void intt(vector<ll>& input) {
		_ntt(input, -1);
		const int n_inv = mod_inv(SZ(input), mod);
		for (auto& x : input) x = x * n_inv % mod;
	}

	vector<ll> convolution(const vector<ll>& a, const vector<ll>& b){
		int ntt_size = 1;
		while (ntt_size < SZ(a) + SZ(b)) ntt_size *= 2;

		vector<ll> _a = a, _b = b;
		_a.resize(ntt_size); _b.resize(ntt_size);

		ntt(_a);
		ntt(_b);

		REP(i, ntt_size) (_a[i] *= _b[i]) %= mod;

		intt(_a);
		return _a;
	}
};

using NTT_1 = NTT<167772161, 3>;
using NTT_2 = NTT<469762049, 3>;
using NTT_3 = NTT<1224736769, 3>;

vector<ll> fast_int32mod_convolution(vector<ll> a, vector<ll> b, int mod){
	for (auto& x : a) x %= mod;
	for (auto& x : b) x %= mod;

	NTT_1 ntt1; NTT_2 ntt2; NTT_3 ntt3;
	assert(ntt1.get_mod() < ntt2.get_mod() && ntt2.get_mod() < ntt3.get_mod());
	auto x = ntt1.convolution(a, b);
	auto y = ntt2.convolution(a, b);
	auto z = ntt3.convolution(a, b);

	const ll m1 = ntt1.get_mod(), m2 = ntt2.get_mod(), m3 = ntt3.get_mod();
	const ll m1_inv_m2 = mod_inv<ll>(m1, m2);
	const ll m12_inv_m3 = mod_inv<ll>(m1 * m2, m3);
	const ll m12_mod = m1 * m2 % mod;
	vector<ll> ret(SZ(x));
	REP(i, SZ(x)){
		ll v1 = (y[i] - x[i]) *  m1_inv_m2 % m2;
		if (v1 < 0) v1 += m2;
		ll v2 = (z[i] - (x[i] + m1 * v1) % m3) * m12_inv_m3 % m3;
		if (v2 < 0) v2 += m3;
		ll constants3 = (x[i] + m1 * v1 + m12_mod * v2) % mod;
		if (constants3 < 0) constants3 += mod;
		ret[i] = constants3;
	}
	return ret;
}

template<class T> class Polynomial {
public:
  std::vector<T> coef;
  Polynomial() : coef(1, 0) {}
  Polynomial(int N) : coef(N + 1, 0) {}
  Polynomial(std::initializer_list<T> a) : coef(a) {}
  Polynomial(std::string s) {
    using Pit = std::string::const_iterator;

    auto next = [&](Pit& it) {
      do { ++it; } while (it != s.cend() && *it == ' ');
    };
    auto num = [&](Pit& it) -> int {
      int res = 0;
      while (it != s.cend() && std::isdigit(*it)) {
        res = res * 10 + (*it - '0');
        next(it);
      }
      return res;
    };
    auto atom = [&](Pit& it) -> Polynomial<T> {
      if (std::isdigit(*it)) {
        int c = num(it);
        return { c }; // c
      }
      else if (*it == 'x') {
        next(it);
        if (it != s.cend() && *it == '^') {
          next(it);
          int t = num(it);
          Polynomial<T> r(t);
          r[t] = 1;
          return r; // x^e
        }
        else {
          return {0, 1}; // x
        }
      }
      return {};
    };
    std::function<Polynomial<T>(Pit&)> expr;
    auto mono = [&](Pit& it) -> Polynomial<T> {
      if (*it == '(') {
        next(it);
        auto r = expr(it);
        next(it);
        if (it != s.cend() && *it == '^') {
          next(it);
          r ^= num(it);
        }
        return r;
      }
      else {
        return atom(it);
      }
    };
    auto prod = [&](Pit& it) -> Polynomial<T> {
      Polynomial<T> r({ 1 });
      r *= mono(it);
      while (it != s.cend()) {
        if (*it == '*') {
          next(it);
          r *= mono(it);
        }
        else if (*it == '(' || *it == 'x' || std::isdigit(*it)) {
          r *= mono(it);
        }
        else break;
      }
      return r;
    };
    expr = [&](Pit& it) -> Polynomial<T> {
      Polynomial<T> r = prod(it);
      while (it != s.cend() && *it != ')') {
        bool neg = false;
        if (*it == '-') neg = true;
        next(it);
        if (neg) r -= prod(it);
        else r += prod(it);
      }
      return r;
    };

    if (!s.empty() && (s.front() == '+' || s.front() == '-')) s = '0' + s;
    Pit it = s.cbegin();
    *this = expr(it);
  }
  Polynomial(const char* s) : Polynomial(std::string(s)) {}
  Polynomial& operator+=(const Polynomial& r) {
    if (SZ(coef) < SZ(r.coef)) coef.resize(SZ(r.coef));
    REP(i,SZ(r.coef)) coef[i] = safe_mod(coef[i]+r.coef[i],MOD);
    return *this;
  }
  Polynomial& operator-=(const Polynomial& r) {
    if (SZ(coef) < SZ(r.coef)) coef.resize(SZ(r.coef));
    REP(i,SZ(r.coef)) coef[i] = safe_mod(coef[i]-r.coef[i],MOD);
    return *this;
  }
  //Polynomial& operator*=(const Polynomial & r) { // O(N^2)
  //	std::vector<T> c(this->SZ(coef) + SZ(r.coef) - 1);
  //	REP(i, this->SZ(coef)) REP(j, SZ(r.coef)) c[i + j] += this->coef[i] * r.coef[j];
  //	coef = c;
  //	return *this;
  //}

  Polynomial& operator*=(const Polynomial& r) { // O(NlogN)
    this->coef = fast_int32mod_convolution(this->coef, r.coef, MOD);
    return *this;
  }
  Polynomial operator*(const Polynomial& r) { Polynomial res(*this); return res *= r; }
  Polynomial& mul_sparse(const std::map<int, T>& r, int sizeLimit = -1) {
    size_t sz = SZ(coef);
    size_t rsz = r.rbegin()->first;
    coef.resize(sz + rsz, 0);
    for (int i = sz + rsz - 1; i >= 0; --i) {
      for (const auto& p : r) if (p.first != 0) {
        if (i >= p.first) {
          coef[i] = safe_mod(coef[i] + safe_mod(p.second,MOD) * safe_mod(coef[i - p.first],MOD), MOD);
        }
      }
    }
    if (sizeLimit > 0 && SZ(coef) > sizeLimit) {
      coef.resize(sizeLimit);
    }
    return *this;
  }
  Polynomial& div_sparse(const std::map<int, T>& r, int sz = -1) {
    if (sz == -1) sz = SZ(coef);
    REP(i,sz) {
       for (const auto& p : r) if (p.first != 0) {
        if (i >= p.first) {
          coef[i] = safe_mod(coef[i] - safe_mod(p.second,MOD) * safe_mod(coef[i - p.first],MOD), MOD);
        }
      }
    }
    return *this;
  }
  Polynomial& operator^=(long long p) {
    Polynomial x(*this);
    *this = Polynomial(0);
    coef[0] = 1;
    while (p) {
      if (p & 1) (*this) *= x;
      x *= x;
      p >>= 1;
    }
    return *this;
  }
  Polynomial operator^(long long p) { Polynomial res(*this); return res ^= p; }
  T& operator[](size_t i) { return coef[i]; }
  std::vector<T> getCoef() const { return coef; }
};

// get [X^N]P(X)/Q(X) in O(k logk logN) where k=deg(Q)
template<class T> T coef_of_rational_polynomial(ll N, Polynomial<T> P, Polynomial<T> Q) {
  for (; N; N >>= 1) {
    Polynomial<T> Q1(Q);
    for (int i = 1; i < SZ(Q.coef); i += 2) Q1[i] = MOD-Q[i];

    P *= Q1;
    int j = 0;
    for (int i = N & 1; i < SZ(P.coef); i += 2) P.coef[j++] = P.coef[i];
    P.coef.resize(j);

    Q *= Q1;
    j = 0;
    for (int i = 0; i < SZ(Q.coef); i += 2) Q.coef[j++] = Q.coef[i];
    Q.coef.resize(j);
  }
  return safe_mod(P[0],MOD);
}


#endif
yukicoder

結果

ソースコード
