fork download
copy 
  1. #include <bits/stdc++.h>
  2.  
  3. using namespace std;
  4.  
  5. typedef long long ll;
  6.  
  7. #ifdef iq
  8.   mt19937 rnd(228);
  9. #else
  10.   mt19937 rnd(chrono::high_resolution_clock::now().time_since_epoch().count());
  11. #endif
  12.  
  13.  
  14. template <typename T>
  15. T inverse(T a, T m) {
  16.   T u = 0, v = 1;
  17.   while (a != 0) {
  18.     T t = m / a;
  19.     m -= t * a; swap(a, m);
  20.     u -= t * v; swap(u, v);
  21.   }
  22.   assert(m == 1);
  23.   return u;
  24. }
  25. template <typename T>
  26. class Modular {
  27.  public:
  28.   using Type = typename decay<decltype(T::value)>::type;
  29.   constexpr Modular() : value() {}
  30.   template <typename U>
  31.   Modular(const U& x) {
  32.     value = normalize(x);
  33.   }
  34.   template <typename U>
  35.   static Type normalize(const U& x) {
  36.     Type v;
  37.     if (-mod() <= x && x < mod()) v = static_cast<Type>(x);
  38.     else v = static_cast<Type>(x % mod());
  39.     if (v < 0) v += mod();
  40.     return v;
  41.   }
  42.   const Type& operator()() const { return value; }
  43.   template <typename U>
  44.   explicit operator U() const { return static_cast<U>(value); }
  45.   constexpr static Type mod() { return T::value; }
  46.   Modular& operator+=(const Modular& other) { if ((value += other.value) >= mod()) value -= mod(); return *this; }
  47.   Modular& operator-=(const Modular& other) { if ((value -= other.value) < 0) value += mod(); return *this; }
  48.   template <typename U> Modular& operator+=(const U& other) { return *this += Modular(other); }
  49.   template <typename U> Modular& operator-=(const U& other) { return *this -= Modular(other); }
  50.   Modular& operator++() { return *this += 1; }
  51.   Modular& operator--() { return *this -= 1; }
  52.   Modular operator++(int) { Modular result(*this); *this += 1; return result; }
  53.   Modular operator--(int) { Modular result(*this); *this -= 1; return result; }
  54.   Modular operator-() const { return Modular(-value); }
  55.   template <typename U = T>
  56.   typename enable_if<is_same<typename Modular<U>::Type, int>::value, Modular>::type& operator*=(const Modular& rhs) {
  57. #ifdef _WIN32
  58.     uint64_t x = static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value);
  59.     uint32_t xh = static_cast<uint32_t>(x >> 32), xl = static_cast<uint32_t>(x), d, m;
  60.     asm(
  61.       "divl %4; \n\t"
  62.       : "=a" (d), "=d" (m)
  63.       : "d" (xh), "a" (xl), "r" (mod())
  64.     );
  65.     value = m;
  66. #else
  67.     value = normalize(static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value));
  68. #endif
  69.     return *this;
  70.   }
  71.   template <typename U = T>
  72.   typename enable_if<is_same<typename Modular<U>::Type, int64_t>::value, Modular>::type& operator*=(const Modular& rhs) {
  73.     int64_t q = static_cast<int64_t>(static_cast<long double>(value) * rhs.value / mod());
  74.     value = normalize(value * rhs.value - q * mod());
  75.     return *this;
  76.   }
  77.   template <typename U = T>
  78.   typename enable_if<!is_integral<typename Modular<U>::Type>::value, Modular>::type& operator*=(const Modular& rhs) {
  79.     value = normalize(value * rhs.value);
  80.     return *this;
  81.   }
  82.   Modular& operator/=(const Modular& other) { return *this *= Modular(inverse(other.value, mod())); }
  83.   template <typename U>
  84.   friend bool operator==(const Modular<U>& lhs, const Modular<U>& rhs);
  85.   template <typename U>
  86.   friend bool operator<(const Modular<U>& lhs, const Modular<U>& rhs);
  87.   template <typename U>
  88.   friend std::istream& operator>>(std::istream& stream, Modular<U>& number);
  89.  private:
  90.   Type value;
  91. };
  92. template <typename T> bool operator==(const Modular<T>& lhs, const Modular<T>& rhs) { return lhs.value == rhs.value; }
  93. template <typename T, typename U> bool operator==(const Modular<T>& lhs, U rhs) { return lhs == Modular<T>(rhs); }
  94. template <typename T, typename U> bool operator==(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) == rhs; }
  95. template <typename T> bool operator!=(const Modular<T>& lhs, const Modular<T>& rhs) { return !(lhs == rhs); }
  96. template <typename T, typename U> bool operator!=(const Modular<T>& lhs, U rhs) { return !(lhs == rhs); }
  97. template <typename T, typename U> bool operator!=(U lhs, const Modular<T>& rhs) { return !(lhs == rhs); }
  98. template <typename T> bool operator<(const Modular<T>& lhs, const Modular<T>& rhs) { return lhs.value < rhs.value; }
  99. template <typename T> Modular<T> operator+(const Modular<T>& lhs, const Modular<T>& rhs) { return Modular<T>(lhs) += rhs; }
  100. template <typename T, typename U> Modular<T> operator+(const Modular<T>& lhs, U rhs) { return Modular<T>(lhs) += rhs; }
  101. template <typename T, typename U> Modular<T> operator+(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) += rhs; }
  102. template <typename T> Modular<T> operator-(const Modular<T>& lhs, const Modular<T>& rhs) { return Modular<T>(lhs) -= rhs; }
  103. template <typename T, typename U> Modular<T> operator-(const Modular<T>& lhs, U rhs) { return Modular<T>(lhs) -= rhs; }
  104. template <typename T, typename U> Modular<T> operator-(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) -= rhs; }
  105. template <typename T> Modular<T> operator*(const Modular<T>& lhs, const Modular<T>& rhs) { return Modular<T>(lhs) *= rhs; }
  106. template <typename T, typename U> Modular<T> operator*(const Modular<T>& lhs, U rhs) { return Modular<T>(lhs) *= rhs; }
  107. template <typename T, typename U> Modular<T> operator*(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) *= rhs; }
  108. template <typename T> Modular<T> operator/(const Modular<T>& lhs, const Modular<T>& rhs) { return Modular<T>(lhs) /= rhs; }
  109. template <typename T, typename U> Modular<T> operator/(const Modular<T>& lhs, U rhs) { return Modular<T>(lhs) /= rhs; }
  110. template <typename T, typename U> Modular<T> operator/(U lhs, const Modular<T>& rhs) { return Modular<T>(lhs) /= rhs; }
  111. template<typename T, typename U>
  112. Modular<T> power(const Modular<T>& a, const U& b) {
  113.   assert(b >= 0);
  114.   Modular<T> x = a, res = 1;
  115.   U p = b;
  116.   while (p > 0) {
  117.     if (p & 1) res *= x;
  118.     x *= x;
  119.     p >>= 1;
  120.   }
  121.   return res;
  122. }
  123. template <typename T>
  124. string to_string(const Modular<T>& number) {
  125.   return to_string(number());
  126. }
  127. template <typename T>
  128. std::ostream& operator<<(std::ostream& stream, const Modular<T>& number) {
  129.   return stream << number();
  130. }
  131. template <typename T>
  132. std::istream& operator>>(std::istream& stream, Modular<T>& number) {
  133.   typename common_type<typename Modular<T>::Type, int64_t>::type x;
  134.   stream >> x;
  135.   number.value = Modular<T>::normalize(x);
  136.   return stream;
  137. }
  138. /*
  139. using ModType = int;
  140. struct VarMod { static ModType value; };
  141. ModType VarMod::value;
  142. ModType& md = VarMod::value;
  143. using Mint = Modular<VarMod>;
  144. */
  145. constexpr int md = 998244353;
  146. using Mint = Modular<std::integral_constant<decay<decltype(md)>::type, md>>;
  147.  
  148.   template <typename T>
  149. class NTT {
  150.  public:
  151.   using Type = typename decay<decltype(T::value)>::type;
  152.   static Type md;
  153.   static Modular<T> root;
  154.   static int base;
  155.   static int max_base;
  156.   static vector<Modular<T>> roots;
  157.   static vector<int> rev;
  158.   static void clear() {
  159.     root = 0;
  160.     base = 0;
  161.     max_base = 0;
  162.     roots.clear();
  163.     rev.clear();
  164.   }
  165.   static void init() {
  166.     md = T::value;
  167.     assert(md >= 3 && md % 2 == 1);
  168.     auto tmp = md - 1;
  169.     max_base = 0;
  170.     while (tmp % 2 == 0) {
  171.       tmp /= 2;
  172.       max_base++;
  173.     }
  174.     root = 2;
  175.     while (power(root, (md - 1) >> 1) == 1) {
  176.       root++;
  177.     }
  178.     assert(power(root, md - 1) == 1);
  179.     root = power(root, (md - 1) >> max_base);
  180.     base = 1;
  181.     rev = {0, 1};
  182.     roots = {0, 1};
  183.   }
  184.   static void ensure_base(int nbase) {
  185.     if (md != T::value) {
  186.       clear();
  187.     }
  188.     if (roots.empty()) {
  189.       init();
  190.     }
  191.     if (nbase <= base) {
  192.       return;
  193.     }
  194.     assert(nbase <= max_base);
  195.     rev.resize(1 << nbase);
  196.     for (int i = 0; i < (1 << nbase); i++) {
  197.       rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (nbase - 1));
  198.     }
  199.     roots.resize(1 << nbase);
  200.     while (base < nbase) {
  201.       Modular<T> z = power(root, 1 << (max_base - 1 - base));
  202.       for (int i = 1 << (base - 1); i < (1 << base); i++) {
  203.         roots[i << 1] = roots[i];
  204.         roots[(i << 1) + 1] = roots[i] * z;
  205.       }
  206.       base++;
  207.     }
  208.   }
  209.   static void fft(vector<Modular<T>> &a) {
  210.     int n = (int) a.size();
  211.     assert((n & (n - 1)) == 0);
  212.     int zeros = __builtin_ctz(n);
  213.     ensure_base(zeros);
  214.     int shift = base - zeros;
  215.     for (int i = 0; i < n; i++) {
  216.       if (i < (rev[i] >> shift)) {
  217.         swap(a[i], a[rev[i] >> shift]);
  218.       }
  219.     }
  220.     for (int k = 1; k < n; k <<= 1) {
  221.       for (int i = 0; i < n; i += 2 * k) {
  222.         for (int j = 0; j < k; j++) {
  223.           Modular<T> x = a[i + j];
  224.           Modular<T> y = a[i + j + k] * roots[j + k];
  225.           a[i + j] = x + y;
  226.           a[i + j + k] = x - y;
  227.         }
  228.       }
  229.     }
  230.   }
  231.   static vector<Modular<T>> multiply(vector<Modular<T>> a, vector<Modular<T>> b) {
  232.     if (a.empty() || b.empty()) {
  233.       return {};
  234.     }
  235.     int eq = (a == b);
  236.     int need = (int) a.size() + (int) b.size() - 1;
  237.     int nbase = 0;
  238.     while ((1 << nbase) < need) nbase++;
  239.     ensure_base(nbase);
  240.     int sz = 1 << nbase;
  241.     a.resize(sz);
  242.     b.resize(sz);
  243.     fft(a);
  244.     if (eq) b = a; else fft(b);
  245.     Modular<T> inv_sz = 1 / static_cast<Modular<T>>(sz);
  246.     for (int i = 0; i < sz; i++) {
  247.       a[i] *= b[i] * inv_sz;
  248.     }
  249.     reverse(a.begin() + 1, a.end());
  250.     fft(a);
  251.     a.resize(need);
  252.     return a;
  253.   }
  254. };
  255. template <typename T> typename NTT<T>::Type NTT<T>::md;
  256. template <typename T> Modular<T> NTT<T>::root;
  257. template <typename T> int NTT<T>::base;
  258. template <typename T> int NTT<T>::max_base;
  259. template <typename T> vector<Modular<T>> NTT<T>::roots;
  260. template <typename T> vector<int> NTT<T>::rev;
  261. template <typename T>
  262. vector<Modular<T>> inverse(const vector<Modular<T>>& a) {
  263.   assert(!a.empty());
  264.   int n = (int) a.size();
  265.   vector<Modular<T>> b = {1 / a[0]};
  266.   while ((int) b.size() < n) {
  267.     vector<Modular<T>> x(a.begin(), a.begin() + min(a.size(), b.size() << 1));
  268.     x.resize(b.size() << 1);
  269.     b.resize(b.size() << 1);
  270.     vector<Modular<T>> c = b;
  271.     NTT<T>::fft(c);
  272.     NTT<T>::fft(x);
  273.     Modular<T> inv = 1 / static_cast<Modular<T>>((int) x.size());
  274.     for (int i = 0; i < (int) x.size(); i++) {
  275.       x[i] *= c[i] * inv;
  276.     }
  277.     reverse(x.begin() + 1, x.end());
  278.     NTT<T>::fft(x);
  279.     rotate(x.begin(), x.begin() + (x.size() >> 1), x.end());
  280.     fill(x.begin() + (x.size() >> 1), x.end(), 0);
  281.     NTT<T>::fft(x);
  282.     for (int i = 0; i < (int) x.size(); i++) {
  283.       x[i] *= c[i] * inv;
  284.     }
  285.     reverse(x.begin() + 1, x.end());
  286.     NTT<T>::fft(x);
  287.     for (int i = 0; i < ((int) x.size() >> 1); i++) {
  288.       b[i + ((int) x.size() >> 1)] = -x[i];
  289.     }
  290.   }
  291.   b.resize(n);
  292.   return b;
  293. }
  294. template <typename T>
  295. vector<Modular<T>> inverse_old(vector<Modular<T>> a) {
  296.   assert(!a.empty());
  297.   int n = (int) a.size();
  298.   if (n == 1) {
  299.     return {1 / a[0]};
  300.   }
  301.   int m = (n + 1) >> 1;
  302.   vector<Modular<T>> b = inverse(vector<Modular<T>>(a.begin(), a.begin() + m));
  303.   int need = n << 1;
  304.   int nbase = 0;
  305.   while ((1 << nbase) < need) {
  306.     ++nbase;
  307.   }
  308.   NTT<T>::ensure_base(nbase);
  309.   int size = 1 << nbase;
  310.   a.resize(size);
  311.   b.resize(size);
  312.   NTT<T>::fft(a);
  313.   NTT<T>::fft(b);
  314.   Modular<T> inv = 1 / static_cast<Modular<T>>(size);
  315.   for (int i = 0; i < size; ++i) {
  316.     a[i] = (2 - a[i] * b[i]) * b[i] * inv;
  317.   }
  318.   reverse(a.begin() + 1, a.end());
  319.   NTT<T>::fft(a);
  320.   a.resize(n);
  321.   return a;
  322. }
  323. template <typename T>
  324. vector<Modular<T>> operator*(const vector<Modular<T>>& a, const vector<Modular<T>>& b) {
  325.   if (a.empty() || b.empty()) {
  326.     return {};
  327.   }
  328.   if (min(a.size(), b.size()) < 150) {
  329.     vector<Modular<T>> c(a.size() + b.size() - 1, 0);
  330.     for (int i = 0; i < (int) a.size(); i++) {
  331.       for (int j = 0; j < (int) b.size(); j++) {
  332.         c[i + j] += a[i] * b[j];
  333.       }
  334.     }
  335.     return c;
  336.   }
  337.   return NTT<T>::multiply(a, b);
  338. }
  339. template <typename T>
  340. vector<Modular<T>>& operator*=(vector<Modular<T>>& a, const vector<Modular<T>>& b) {
  341.   return a = a * b;
  342. }
  343.  
  344. const int N = 1e6 + 228;
  345.  
  346. Mint fact[N], rfact[N];
  347.  
  348. int main() {
  349. #ifdef iq
  350.   freopen("a.in", "r", stdin);
  351. #endif
  352.   ios::sync_with_stdio(0);
  353.   cin.tie(0);
  354.   fact[0] = 1;
  355.   for (int i = 1; i < N; i++) {
  356.     fact[i] = fact[i - 1] * i;
  357.   }
  358.   rfact[N - 1] = 1 / fact[N - 1];
  359.   for (int i = N - 2; i >= 0; i--) {
  360.     rfact[i] = (i + 1) * rfact[i + 1];
  361.   }
  362.   int t;
  363.   cin >> t;
  364.   while (t--) {
  365.     int n, p, r;
  366.     cin >> n >> p >> r;
  367.     //x(x-1)...(x-n+1)
  368.     function<vector<Mint>(int)> f = [&] (int n) {
  369.       if (n == 0) {
  370.         vector <Mint> t;
  371.         t.push_back(1);
  372.         return t;
  373.       } else if (n % 2) {
  374.         auto ok = f(n - 1);
  375.         vector <Mint> p;
  376.         p.push_back(-n + 1);
  377.         p.push_back(1);
  378.         return ok * p;
  379.       } else {
  380.         auto ok = f(n / 2);
  381.         int c = -(n / 2);
  382.         //((x+c)^i) -> sum(x^j*c^(i-j)*i!/j!/(i-j)!)
  383.         auto a = ok;
  384.         for (int i = 0; i < (int) a.size(); i++) {
  385.           a[i] = a[i] * fact[i];
  386.         }
  387.         vector <Mint> t(ok.size());
  388.         t[0] = 1;
  389.         for (int i = 1; i < (int) ok.size(); i++) {
  390.           t[i] = t[i - 1] * c;
  391.         }
  392.         for (int i = 0; i < (int) t.size(); i++) {
  393.           t[i] = t[i] * rfact[i];
  394.         }
  395.         reverse(t.begin(), t.end());
  396.         auto mda = a * t;
  397.         int sign = (int) t.size() - 1;
  398.         vector <Mint> b(mda.size() - sign);
  399.         for (int i = sign; i < (int) mda.size(); i++) {
  400.           b[i - sign] += mda[i];
  401.         }
  402.         for (int i = 0; i < (int) b.size(); i++) {
  403.           b[i] = b[i] * rfact[i];
  404.         }
  405.         return ok * b;
  406.       }
  407.     };
  408.     auto ans = f(r);
  409.     n++;
  410.     Mint me = 1;
  411.     Mint ret = 0;
  412.     for (int i = 0; i < (int) ans.size(); i++) {
  413.       if (me == 1) {
  414.         ret += n * ans[i];
  415.       } else {
  416.         ret += (power(me, n) - 1) / (me - 1) * ans[i];
  417.       }
  418.       /*
  419.       Ot cur = power(me, n) - 1;
  420.       // += ans[i] * (me^n - 1)/(me-1)
  421.       //me^0 + me^1 + ... +me^(n-1)
  422.       */
  423.       me *= p;
  424.     }
  425.     cout << ret * rfact[r] << endl;
  426.   }
  427. }
Success #stdin #stdout 0.02s 11380KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/FUI92w
language:
C++14 (gcc 8.3)
created:
4 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!