API tools faq
paste
Advertisement
EntropyIncreaser

Polynomial & Number Theory

EntropyIncreaser
Dec 19th, 2018
1,364
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 18.78 KB | None | 0 0
raw download clone embed print report
  1. #include <cstdio>
  2. #include <cstring>
  3.  
  4. #include <algorithm>
  5. #include <chrono>
  6. #include <random>
  7. #include <functional>
  8. #include <vector>
  9.  
  10. #define LOG(FMT...) fprintf(stderr, FMT)
  11.  
  12. using namespace std;
  13.  
  14. typedef long long ll;
  15.  
  16. const int P = 998244353, R = 3;
  17. const int BRUTE_N2_LIMIT = 50;
  18.  
  19. int mpow(int x, int k, int p = P) {
  20.   int ret = 1;
  21.   while (k) {
  22.     if (k & 1)
  23.       ret = ret * (ll)x % p;
  24.     x = x * (ll)x % p;
  25.     k >>= 1;
  26.   }
  27.   return ret;
  28. }
  29.  
  30. int norm(int x) { return x >= P ? x - P : x; }
  31.  
  32. struct NumberTheory {
  33.  
  34.   typedef pair<int, int> _P2_Field;
  35.  
  36.   mt19937 rng;
  37.  
  38.   NumberTheory() : rng(chrono::steady_clock::now().time_since_epoch().count()) {}
  39.  
  40.   void _exGcd(int a, int b, int& x, int& y) {
  41.     if (!b) {
  42.       x = 1;
  43.       y = 0;
  44.       return;
  45.     }
  46.     _exGcd(b, a % b, y, x);
  47.     y -= a / b * x;
  48.   }
  49.  
  50.   int inv(int a, int p = P) {
  51.     int x, y;
  52.     _exGcd(a, p, x, y);
  53.     if (x < 0)
  54.       x += p;
  55.     return x;
  56.   }
  57.  
  58.   template <class Integer>
  59.   bool quadRes(Integer a, Integer b) {
  60.     if (a <= 1)
  61.       return true;
  62.     while (a % 4 == 0)
  63.       a /= 4;
  64.     if (a % 2 == 0)
  65.       return (b % 8 == 1 || b % 8 == 7) == quadRes(a / 2, b);
  66.     return ((a - 1) % 4 == 0 || (b - 1) % 4 == 0) == quadRes(b % a, a);
  67.   }
  68.  
  69.   // assume p in prime, x in quadratic residue
  70.   int sqrt(int x, int p = P) {
  71.     if (p == 2 || x <= 1)
  72.       return x;
  73.     int w, v, k = (p + 1) / 2;
  74.     do {
  75.       w = rng() % p;
  76.     } while (quadRes(v = int((w * (ll)w - x + p) % p), p));
  77.     _P2_Field res(1, 0), a(w, 1);
  78.     while (k) {
  79.       if (k & 1)
  80.         res = _P2_Field((res.first * (ll)a.first + res.second * (ll)a.second % p * v) % p, (res.first * (ll)a.second + res.second * (ll)a.first) % p);
  81.       if (k >>= 1)
  82.         a = _P2_Field((a.first * (ll)a.first + a.second * (ll)a.second % p * v) % p, (a.first * (ll)a.second << 1) % p);
  83.     }
  84.     return min(res.first, p - res.first);
  85.   }
  86.  
  87. } nt;
  88.  
  89. template <class T, class Comp>
  90. struct AdditionChain {
  91.   int k;
  92.   vector<T> prepare;
  93.   T t, unit;
  94.   Comp comp;
  95.  
  96.   AdditionChain(const T& t, const Comp& comp, int k, const T& unit = 1) : comp(comp), t(t), unit(unit), k(k), prepare(1U << k) {
  97.     prepare[0] = unit;
  98.     for (int i = 1; i < 1 << k; ++i)
  99.       prepare[i] = comp(prepare[i - 1], t);
  100.   }
  101.  
  102.   static AdditionChain fourRussians(const T &t, const Comp &comp, int lgn, const T &unit = 1) {
  103.     lgn = max(lgn, 1);
  104.     int k = 1, lglgn = 1;
  105.     while (2 << lglgn <= lgn)
  106.       ++lglgn;
  107.     int w = lglgn / lgn;
  108.     while (1 << k < w)
  109.       ++k;
  110.     return AdditionChain(t, comp, k, unit);
  111.   }
  112.  
  113.   T pow(int n) const {
  114.     if (n < 1 << k)
  115.       return prepare[n];
  116.     int r = n & ((1 << k) - 1);
  117.     T step = pow(n >> k);
  118.     for (int rep = 0; rep < k; ++rep)
  119.       step = comp(step, step);
  120.     return comp(step, prepare[r]);
  121.   }
  122. };
  123.  
  124. struct Simple {
  125.   int n;
  126.   vector<int> fac, ifac, inv;
  127.  
  128.   void build(int n) {
  129.     this->n = n;
  130.     fac.resize(n + 1);
  131.     ifac.resize(n + 1);
  132.     inv.resize(n + 1);
  133.     fac[0] = 1;
  134.     for (int x = 1; x <= n; ++x)
  135.       fac[x] = fac[x - 1] * (ll)x % P;
  136.     inv[1] = 1;
  137.     for (int x = 2; x <= n; ++x)
  138.       inv[x] = -(P / x) * (ll)inv[P % x] % P + P;
  139.     ifac[0] = 1;
  140.     for (int x = 1; x <= n; ++x)
  141.       ifac[x] = ifac[x - 1] * (ll)inv[x] % P;
  142.   }
  143.  
  144.   Simple() {
  145.     build(1);
  146.   }
  147.  
  148.   void check(int k) {
  149.     int nn = n;
  150.     if (k > nn) {
  151.       while (k > nn)
  152.         nn <<= 1;
  153.       build(nn);
  154.     }
  155.   }
  156.  
  157.   int gfac(int k) {
  158.     check(k);
  159.     return fac[k];
  160.   }
  161.  
  162.   int gifac(int k) {
  163.     check(k);
  164.     return ifac[k];
  165.   }
  166.  
  167.   int ginv(int k) {
  168.     check(k);
  169.     return inv[k];
  170.   }
  171.  
  172.   int binom(int n, int m) {
  173.     if (m < 0 || m > n)
  174.       return 0;
  175.     return gfac(n) * (ll)gifac(m) % P * gifac(n - m) % P;
  176.   }
  177. } simp;
  178.  
  179. const int L2 = 11;
  180.  
  181. struct NTT {
  182.   int L;
  183.   int brev[1 << L2];
  184.   vector<int> root;
  185.  
  186.   NTT() : L(-1) {
  187.     for (int i = 1; i < (1 << L2); ++i)
  188.       brev[i] = brev[i >> 1] >> 1 | ((i & 1) << (L2 - 1));
  189.   }
  190.  
  191.   void prepRoot(int l) {
  192.     L = l;
  193.     root.resize(1 << L);
  194.     int n = 1 << L;
  195.     int primitive = mpow(R, (P - 1) >> L);
  196.     root[0] = 1;
  197.     for (int i = 1; i < n; ++i) root[i] = root[i - 1] * (ll)primitive % P;
  198.   }
  199.  
  200.   void fft(int* a, int lgn, int d = 1) {
  201.     if (L < lgn) prepRoot(lgn);
  202.     int n = 1 << lgn;
  203.     for (int i = 0; i < n; ++i) {
  204.       int rev = (brev[i >> L2] | (brev[i & ((1 << L2) - 1)] << L2)) >> ((L2 << 1) - lgn);
  205.       if (i < rev)
  206.         swap(a[i], a[rev]);
  207.     }
  208.     int rt = d == 1 ? R : nt.inv(R);
  209.     for (int k = L - 1, t = 1; t < n; t <<= 1, --k) {
  210.       for (int i = 0; i < n; i += t << 1) {
  211.         int *p1 = a + i, *p2 = a + i + t;
  212.         for (int j = 0; j < t; ++j) {
  213.           int x = p2[j] * (ll)root[j << k] % P;
  214.           p2[j] = norm(p1[j] + P - x);
  215.           p1[j] = norm(p1[j] + x);
  216.         }
  217.       }
  218.     }
  219.     if (d == -1) {
  220.       reverse(a + 1, a + n);
  221.       int nv = mpow(n, P - 2);
  222.       for (int i = 0; i < n; ++i) a[i] = a[i] * (ll)nv % P;
  223.     }
  224.   }
  225. } ntt;
  226.  
  227. struct Poly {
  228.   vector<int> a;
  229.  
  230.   Poly(int v = 0) : a(1) {
  231.     if ((v %= P) < 0)
  232.       v += P;
  233.     a[0] = v;
  234.   }
  235.  
  236.   Poly(const vector<int>& a) : a(a) {}
  237.  
  238.   Poly(initializer_list<int> init) : a(init) {}
  239.  
  240.   // Helps
  241.   int operator[](int k) const { return k < a.size() ? a[k] : 0; }
  242.  
  243.   int& operator[](int k) {
  244.     if (k >= a.size())
  245.       a.resize(k + 1);
  246.     return a[k];
  247.   }
  248.  
  249.   int deg() const { return a.size() - 1; }
  250.  
  251.   void redeg(int d) { a.resize(d + 1); }
  252.  
  253.   Poly monic() const;
  254.   Poly sunic() const;
  255.  
  256.   Poly slice(int d) const {
  257.     if (d < a.size())
  258.       return vector<int>(a.begin(), a.begin() + d + 1);
  259.     vector<int> res(a);
  260.     res.resize(d + 1);
  261.     return res;
  262.   }
  263.  
  264.   int* base() { return a.begin().base(); }
  265.   const int* base() const { return a.begin().base(); }
  266.  
  267.   Poly println(FILE* fp) const {
  268.     fprintf(fp, "%d", a[0]);
  269.     for (int i = 1; i < a.size(); ++i)
  270.       fprintf(fp, " %d", a[i]);
  271.     fputc('\n', fp);
  272.     return *this;
  273.   }
  274.  
  275.   // Calculations
  276.   Poly operator+(const Poly& rhs) const {
  277.     vector<int> res(max(a.size(), rhs.a.size()));
  278.     for (int i = 0; i < res.size(); ++i)
  279.       if ((res[i] = operator[](i) + rhs[i]) >= P)
  280.         res[i] -= P;
  281.     return res;
  282.   }
  283.  
  284.   Poly operator-() const {
  285.     Poly ret(a);
  286.     for (int i = 0; i < a.size(); ++i)
  287.       if (ret[i])
  288.         ret[i] = P - ret[i];
  289.     return ret;
  290.   }
  291.  
  292.   Poly operator-(const Poly& rhs) const { return operator+(-rhs); }
  293.  
  294.   Poly operator*(const Poly& rhs) const;
  295.  
  296.   Poly operator/(const Poly& rhs) const;
  297.  
  298.   Poly operator%(const Poly& rhs) const;
  299.  
  300.   Poly der() const; // default: remove trailing
  301.   Poly integ() const;
  302.  
  303.   Poly inv() const;
  304.   Poly sqrt() const;
  305.   Poly ln() const;
  306.   Poly exp() const;
  307.   pair<Poly, Poly> sqrti() const;
  308.   pair<Poly, Poly> expi() const;
  309.  
  310.   Poly quo(const Poly& rhs) const;
  311.   pair<Poly, Poly> iquo(const Poly& rhs) const;
  312.  
  313.   pair<Poly, Poly> div(const Poly& rhs) const;
  314.  
  315.   Poly taylor(int k) const;
  316.   Poly pow(int k) const;
  317.   Poly exp(int k) const;
  318. };
  319.  
  320. Poly zeroes(int deg) { return vector<int>(deg + 1); }
  321.  
  322. struct Newton {
  323.   void inv(const Poly& f, const Poly& nttf, Poly& g, const Poly& nttg, int t) {
  324.     int n = 1 << t;
  325.     Poly prod = nttf;
  326.     for (int i = 0; i < (n << 1); ++i)
  327.       prod[i] = prod[i] * (ll)nttg[i] % P;
  328.     ntt.fft(prod.base(), t + 1, -1);
  329.     for (int i = 0; i < n; ++i)
  330.       prod[i] = 0;
  331.     ntt.fft(prod.base(), t + 1, 1);
  332.     for (int i = 0; i < (n << 1); ++i)
  333.       prod[i] = prod[i] * (ll)nttg[i] % P;
  334.     ntt.fft(prod.base(), t + 1, -1);
  335.     for (int i = 0; i < n; ++i)
  336.       prod[i] = 0;
  337.     g = g - prod;
  338.   }
  339.  
  340.   void inv(const Poly& f, const Poly& nttf, Poly& g, int t) {
  341.     Poly nttg = g;
  342.     nttg.redeg((2 << t) - 1);
  343.     ntt.fft(nttg.base(), t + 1, 1);
  344.     inv(f, nttf, g, nttg, t);
  345.   }
  346.  
  347.   void inv(const Poly& f, Poly& g, int t) {
  348.     Poly nttg = g;
  349.     nttg.redeg((2 << t) - 1);
  350.     ntt.fft(nttg.base(), t + 1, 1);
  351.     Poly nttf = f;
  352.     nttf.redeg((2 << t) - 1);
  353.     ntt.fft(nttf.base(), t + 1, 1);
  354.     inv(f, nttf, g, nttg, t);
  355.   }
  356.  
  357.   void sqrt(const Poly& f, Poly& g, Poly& nttg, Poly& h, int t) {
  358.     for (int i = 0; i < (1 << t); ++i)
  359.       nttg[i] = mpow(nttg[i], 2);
  360.     ntt.fft(nttg.base(), t, -1);
  361.     nttg = nttg - f;
  362.     for (int i = 0; i < (1 << t); ++i)
  363.       if ((nttg[i + (1 << t)] += nttg[i]) >= P)
  364.         nttg[i + (1 << t)] -= P;
  365.     memset(nttg.base(), 0, sizeof(int) << t);
  366.     ntt.fft(nttg.base(), t + 1, 1);
  367.     Poly tmp = h;
  368.     tmp.redeg((2 << t) - 1);
  369.     ntt.fft(tmp.base(), t + 1, 1);
  370.     for (int i = 0; i < (2 << t); ++i)
  371.       tmp[i] = tmp[i] * (ll)nttg[i] % P;
  372.     ntt.fft(tmp.base(), t + 1, -1);
  373.     memset(tmp.base(), 0, sizeof(int) << t);
  374.     g = g - tmp * nt.inv(2);
  375.   }
  376.  
  377.   void exp(const Poly& f, Poly& g, Poly& nttg, Poly& h, int t) {
  378.     Poly ntth(h);
  379.     ntt.fft(ntth.base(), t, 1);
  380.     Poly dg = g.der().slice((1 << t) - 1);
  381.     ntt.fft(dg.base(), t, 1);
  382.     Poly tmp = zeroes((1 << t) - 1);
  383.     for (int i = 0; i < (1 << t); ++i) {
  384.       tmp[i] = nttg[i << 1] * (ll)ntth[i] % P;
  385.       dg[i] = dg[i] * (ll)ntth[i] % P;
  386.     }
  387.     ntt.fft(tmp.base(), t, -1);
  388.     ntt.fft(dg.base(), t, -1);
  389.     if (--tmp[0] < 0)
  390.       tmp[0] = P - 1;
  391.     dg.redeg((2 << t) - 1);
  392.     Poly df0 = f.der().slice((1 << t) - 1);
  393.     df0[(1 << t) - 1] = 0;
  394.     for (int i = 0; i < (1 << t); ++i) {
  395.       if ((dg[i | 1 << t] = dg[i] - df0[i]) < 0)
  396.         dg[i | 1 << t] += P;
  397.     }
  398.     memcpy(dg.base(), df0.base(), sizeof(int) * ((1 << t) - 1));
  399.     tmp.redeg((2 << t) - 1);
  400.     ntt.fft(tmp.base(), t + 1, 1);
  401.     df0.redeg((2 << t) - 1);
  402.     ntt.fft(df0.base(), t + 1, 1);
  403.     for (int i = 0; i < (2 << t); ++i)
  404.       df0[i] = df0[i] * (ll)tmp[i] % P;
  405.     ntt.fft(df0.base(), t + 1, -1);
  406.     memcpy(df0.base() + (1 << t), df0.base(), sizeof(int) << t);
  407.     memset(df0.base(), 0, sizeof(int) << t);
  408.     dg = (dg - df0).integ().slice((2 << t) - 1) - f;
  409.     ntt.fft(dg.base(), t + 1, 1);
  410.     for (int i = 0; i < (2 << t); ++i)
  411.       tmp[i] = dg[i] * (ll)nttg[i] % P;
  412.     ntt.fft(tmp.base(), t + 1, -1);
  413.     g.redeg((2 << t) - 1);
  414.     for (int i = 1 << t; i < (2 << t); ++i)
  415.       if (tmp[i])
  416.         g[i] = P - tmp[i];
  417.   }
  418. } nit;
  419.  
  420. struct Transposition {
  421.  
  422.   vector<int> _mul(int l, vector<int> res, const Poly& b) {
  423.     vector<int> tmp(1 << l);
  424.     memcpy(tmp.begin().base(), b.a.begin().base(), sizeof(int) * (b.deg() + 1));
  425.     reverse(tmp.begin() + 1, tmp.end());
  426.     ntt.fft(tmp.begin().base(), l, 1);
  427.     for (int i = 0; i < (1 << l); ++i)
  428.       res[i] = res[i] * (ll)tmp[i] % P;
  429.     ntt.fft(res.begin().base(), l, -1);
  430.     return res;
  431.   }
  432.  
  433.   Poly bfMul(const Poly& a, const Poly& b) {
  434.     int n = a.deg(), m = b.deg();
  435.     Poly ret = zeroes(n - m);
  436.     for (int i = 0; i <= n - m; ++i)
  437.       for (int j = 0; j <= m; ++j)
  438.         ret[i] = (ret[i] + a[i + j] * (ll)b[j]) % P;
  439.     return ret;
  440.   }
  441.  
  442.   Poly mul(const Poly& a, const Poly& b) {
  443.     if (a.deg() < b.deg()) return 0;
  444.     if (a.deg() <= BRUTE_N2_LIMIT) return bfMul(a, b);
  445.     int l = 0;
  446.     while ((1 << l) <= a.deg()) ++l;
  447.     vector<int> res(1 << l);
  448.     memcpy(res.begin().base(), a.a.begin().base(), sizeof(int) * (a.deg() + 1));
  449.     ntt.fft(res.begin().base(), l, 1);
  450.     res = _mul(l, res, b);
  451.     res.resize(a.deg() - b.deg() + 1);
  452.     return res;
  453.   }
  454.  
  455.   pair<Poly, Poly> mul2(const Poly& a, const Poly& b1, const Poly& b2) {
  456.     if (a.deg() <= BRUTE_N2_LIMIT) return make_pair(bfMul(a, b1), bfMul(a, b2));
  457.     int l = 0;
  458.     while ((1 << l) <= a.deg()) ++l;
  459.     vector<int> fa(1 << l);
  460.     memcpy(fa.begin().base(), a.a.begin().base(), sizeof(int) * (a.deg() + 1));
  461.     ntt.fft(fa.begin().base(), l, 1);
  462.     vector<int> res1 = _mul(l, fa, b1), res2 = _mul(l, fa, b2);
  463.     res1.resize(a.deg() - b1.deg() + 1);
  464.     res2.resize(a.deg() - b2.deg() + 1);
  465.     return make_pair(res1, res2);
  466.   }
  467.  
  468.   vector<int> ls, rs, pos;
  469.   vector<Poly> p, q;
  470.  
  471.   void _build(int n) {
  472.     ls.assign(n * 2 - 1, 0);
  473.     rs.assign(n * 2 - 1, 0);
  474.     p.assign(n * 2 - 1, 0);
  475.     q.assign(n * 2 - 1, 0);
  476.     pos.resize(n);
  477.     int cnt = 0;
  478.     function<int(int, int)> dfs = [&](int l, int r) {
  479.       if (l == r) {
  480.         pos[l] = cnt;
  481.         return cnt++;
  482.       }
  483.       int ret = cnt++;
  484.       int mid = (l + r) >> 1;
  485.       ls[ret] = dfs(l, mid);
  486.       rs[ret] = dfs(mid + 1, r);
  487.       return ret;
  488.     };
  489.     dfs(0, n - 1);
  490.   }
  491.  
  492.   vector<int> _eval(vector<int> f, const vector<int>& x) {
  493.     int n = f.size();
  494.     _build(n);
  495.     for (int i = 0; i < n; ++i)
  496.       q[pos[i]] = {1, norm(P - x[i])};
  497.     for (int i = n * 2 - 2; i >= 0; --i)
  498.       if (ls[i])
  499.         q[i] = q[ls[i]] * q[rs[i]];
  500.     f.resize(n * 2);
  501.     p[0] = mul(f, q[0].inv());
  502.     for (int i = 0; i < n * 2 - 1; ++i)
  503.       if (ls[i])
  504.         tie(p[ls[i]], p[rs[i]]) = mul2(p[i], q[rs[i]], q[ls[i]]);
  505.     vector<int> ret(n);
  506.     for (int i = 0; i < n; ++i)
  507.       ret[i] = p[pos[i]][0];
  508.     return ret;
  509.   }
  510.  
  511.   vector<int> eval(const Poly& f, const vector<int>& x) {
  512.     int n = f.deg() + 1, m = x.size();
  513.     vector<int> tmpf = f.a, tmpx = x;
  514.     tmpf.resize(max(n, m));
  515.     tmpx.resize(max(n, m));
  516.     vector<int> ret = _eval(tmpf, tmpx);
  517.     ret.resize(m);
  518.     return ret;
  519.   }
  520.  
  521.   Poly inter(const vector<int>& x, const vector<int>& y) {
  522.     int n = x.size();
  523.     _build(n);
  524.     for (int i = 0; i < n; ++i)
  525.       q[pos[i]] = {1, norm(P - x[i])};
  526.     for (int i = n * 2 - 2; i >= 0; --i)
  527.       if (ls[i])
  528.         q[i] = q[ls[i]] * q[rs[i]];
  529.     Poly tmp = q[0];
  530.     reverse(tmp.a.begin(), tmp.a.end());
  531.     vector<int> f = tmp.der().a;
  532.     f.resize(n * 2);
  533.     p[0] = mul(f, q[0].inv());
  534.     for (int i = 0; i < n * 2 - 1; ++i)
  535.       if (ls[i])
  536.         tie(p[ls[i]], p[rs[i]]) = mul2(p[i], q[rs[i]], q[ls[i]]);
  537.     for (int i = 0; i < n; ++i)
  538.       p[pos[i]] = nt.inv(p[pos[i]][0]) * (ll)y[i] % P;
  539.     for (int i = 0; i < n * 2 - 1; ++i)
  540.       reverse(q[i].a.begin(), q[i].a.end());
  541.     for (int i = n * 2 - 2; i >= 0; --i)
  542.       if (ls[i])
  543.         p[i] = p[ls[i]] * q[rs[i]] + p[rs[i]] * q[ls[i]];
  544.     return p[0];
  545.   }
  546.  
  547. } tp;
  548.  
  549. Poly operator "" _z(unsigned long long a) { return {0, (int)a}; }
  550.  
  551. Poly operator+(int v, const Poly& rhs) { return Poly(v) + rhs; }
  552.  
  553. Poly Poly::operator*(const Poly& rhs) const {
  554.   int n = deg(), m = rhs.deg();
  555.   if (n <= 10 || m <= 10 || n + m <= BRUTE_N2_LIMIT) {
  556.     Poly ret = zeroes(n + m);
  557.     for (int i = 0; i <= n; ++i)
  558.       for (int j = 0; j <= m; ++j)
  559.         ret[i + j] = (ret[i + j] + a[i] * (ll)rhs[j]) % P;
  560.     return ret;
  561.   }
  562.   n += m;
  563.   int l = 0;
  564.   while ((1 << l) <= n)
  565.     ++l;
  566.   vector<int> res(1 << l), tmp(1 << l);
  567.   memcpy(res.begin().base(), base(), a.size() * sizeof(int));
  568.   ntt.fft(res.begin().base(), l, 1);
  569.   memcpy(tmp.begin().base(), rhs.base(), rhs.a.size() * sizeof(int));
  570.   ntt.fft(tmp.begin().base(), l, 1);
  571.   for (int i = 0; i < (1 << l); ++i)
  572.     res[i] = res[i] * (ll)tmp[i] % P;
  573.   ntt.fft(res.begin().base(), l, -1);
  574.   res.resize(n + 1);
  575.   return res;
  576. }
  577.  
  578. Poly Poly::inv() const {
  579.   Poly g = nt.inv(a[0]);
  580.   for (int t = 0; (1 << t) <= deg(); ++t)
  581.     nit.inv(slice((2 << t) - 1), g, t);
  582.   g.redeg(deg());
  583.   return g;
  584. }
  585.  
  586. Poly Poly::taylor(int k) const {
  587.   int n = deg();
  588.   Poly t = zeroes(n);
  589.   simp.check(n);
  590.   for (int i = 0; i <= n; ++i)
  591.     t[n - i] = a[i] * (ll)simp.fac[i] % P;
  592.   int pw = 1;
  593.   Poly help = vector<int>(simp.ifac.begin(), simp.ifac.begin() + n + 1);
  594.   for (int i = 0; i <= n; ++i) {
  595.     help[i] = help[i] * (ll)pw % P;
  596.     pw = pw * (ll)k % P;
  597.   }
  598.   t = t * help;
  599.   for (int i = 0; i <= n; ++i)
  600.     help[i] = t[n - i] * (ll)simp.ifac[i] % P;
  601.   return help;
  602. }
  603.  
  604. Poly Poly::pow(int k) const {
  605.   if (k == 0)
  606.     return 1;
  607.   if (k == 1)
  608.     return *this;
  609.   int n = deg() * k;
  610.   int lgn = 0;
  611.   while ((1 << lgn) <= n)
  612.     ++lgn;
  613.   vector<int> val = a;
  614.   val.resize(1 << lgn);
  615.   ntt.fft(val.begin().base(), lgn, 1);
  616.   for (int i = 0; i < (1 << lgn); ++i)
  617.     val[i] = mpow(val[i], k);
  618.   ntt.fft(val.begin().base(), lgn, -1);
  619.   return val;
  620. }
  621.  
  622. Poly Poly::der() const {
  623.   if (deg() == 0)
  624.     return 0;
  625.   vector<int> res(deg());
  626.   for (int i = 0; i < deg(); ++i)
  627.     res[i] = a[i + 1] * (ll)(i + 1) % P;
  628.   return res;
  629. }
  630.  
  631. Poly Poly::integ() const {
  632.   vector<int> res(deg() + 2);
  633.   simp.check(deg() + 1);
  634.   for (int i = 0; i <= deg(); ++i)
  635.     res[i + 1] = a[i] * (ll)simp.inv[i + 1] % P;
  636.   return res;
  637. }
  638.  
  639. Poly Poly::quo(const Poly &rhs) const {
  640.   if (rhs.deg() == 0)
  641.     return a[0] * (ll) nt.inv(rhs[0]) % P;
  642.   Poly g = nt.inv(rhs[0]);
  643.   int t = 0, n;
  644.   for (n = 1; (n << 1) <= rhs.deg(); ++t, n <<= 1)
  645.     nit.inv(rhs.slice((n << 1) - 1), g, t);
  646.   Poly nttg = g;
  647.   nttg.redeg((n << 1) - 1);
  648.   ntt.fft(nttg.base(), t + 1, 1);
  649.   Poly eps1 = rhs.slice((n << 1) - 1);
  650.   ntt.fft(eps1.base(), t + 1, 1);
  651.   for (int i = 0; i < (n << 1); ++i)
  652.     eps1[i] = eps1[i] * (ll) nttg[i] % P;
  653.   ntt.fft(eps1.base(), t + 1, -1);
  654.   memcpy(eps1.base(), eps1.base() + n, sizeof(int) << t);
  655.   memset(eps1.base() + n, 0, sizeof(int) << t);
  656.   ntt.fft(eps1.base(), t + 1, 1);
  657.   Poly h0 = slice(n - 1);
  658.   h0.redeg((n << 1) - 1);
  659.   ntt.fft(h0.base(), t + 1);
  660.   Poly h0g0 = zeroes((n << 1) - 1);
  661.   for (int i = 0; i < (n << 1); ++i)
  662.     h0g0[i] = h0[i] * (ll)nttg[i] % P;
  663.   ntt.fft(h0g0.base(), t + 1, -1);
  664.   Poly h0eps1 = zeroes((n << 1) - 1);
  665.   for (int i = 0; i < (n << 1); ++i)
  666.     h0eps1[i] = h0[i] * (ll)eps1[i] % P;
  667.   ntt.fft(h0eps1.base(), t + 1, -1);
  668.   for (int i = 0; i < n; ++i) {
  669.     h0eps1[i] = operator[](i + n) - h0eps1[i];
  670.     if (h0eps1[i] < 0)
  671.       h0eps1[i] += P;
  672.   }
  673.   memset(h0eps1.base() + n, 0, sizeof(int) << t);
  674.   ntt.fft(h0eps1.base(), t + 1);
  675.   for (int i = 0; i < (n << 1); ++i)
  676.     h0eps1[i] = h0eps1[i] * (ll)nttg[i] % P;
  677.   ntt.fft(h0eps1.base(), t + 1, -1);
  678.   memcpy(h0eps1.base() + n, h0eps1.base(), sizeof(int) << t);
  679.   memset(h0eps1.base(), 0, sizeof(int) << t);
  680.   return (h0g0 + h0eps1).slice(rhs.deg());
  681. }
  682.  
  683. Poly Poly::ln() const {
  684.   if (deg() == 0)
  685.     return 0;
  686.   return der().quo(slice(deg() - 1)).integ();
  687. }
  688.  
  689. pair<Poly, Poly> Poly::sqrti() const {
  690.   Poly g = nt.sqrt(a[0]), h = nt.inv(g[0]), nttg = g;
  691.   for (int t = 0; (1 << t) <= deg(); ++t) {
  692.     nit.sqrt(slice((2 << t) - 1), g, nttg, h, t);
  693.     nttg = g;
  694.     ntt.fft(nttg.base(), t + 1, 1);
  695.     nit.inv(g, nttg, h, t);
  696.   }
  697.   return make_pair(g.slice(deg()), h.slice(deg()));
  698. }
  699.  
  700. Poly Poly::sqrt() const {
  701.   Poly g = nt.sqrt(a[0]), h = nt.inv(g[0]), nttg = g;
  702.   for (int t = 0; (1 << t) <= deg(); ++t) {
  703.     nit.sqrt(slice((2 << t) - 1), g, nttg, h, t);
  704.     if ((2 << t) <= deg()) {
  705.       nttg = g;
  706.       ntt.fft(nttg.base(), t + 1, 1);
  707.       nit.inv(g, nttg, h, t);
  708.     }
  709.   }
  710.   return g.slice(deg());
  711. }
  712.  
  713. Poly Poly::exp() const {
  714.   Poly g = 1, h = 1, nttg = {1, 1};
  715.   for (int t = 0; (1 << t) <= deg(); ++t) {
  716.     nit.exp(slice((2 << t) - 1), g, nttg, h, t);
  717.     if ((2 << t) <= deg()) {
  718.       nttg = g;
  719.       nttg.redeg((4 << t) - 1);
  720.       ntt.fft(nttg.base(), t + 2);
  721.       Poly f2n = zeroes((2 << t) - 1);
  722.       for (int i = 0; i < (2 << t); ++i)
  723.         f2n[i] = nttg[i << 1];
  724.       nit.inv(g, f2n, h, t);
  725.     } else {
  726.       nttg = g;
  727.       ntt.fft(nttg.base(), t + 1, 1);
  728.     }
  729.   }
  730.   return g.slice(deg());
  731. }
  732.  
  733. pair<Poly, Poly> Poly::expi() const {
  734.   Poly g = 1, h = 1, nttg = {1, 1};
  735.   for (int t = 0; (1 << t) <= deg(); ++t) {
  736.     nit.exp(slice((2 << t) - 1), g, nttg, h, t);
  737.     nttg = g;
  738.     nttg.redeg((4 << t) - 1);
  739.     ntt.fft(nttg.base(), t + 2);
  740.     Poly f2n = zeroes((2 << t) - 1);
  741.     for (int i = 0; i < (2 << t); ++i)
  742.       f2n[i] = nttg[i << 1];
  743.     nit.inv(g, f2n, h, t);
  744.   }
  745.   return make_pair(g.slice(deg()), h.slice(deg()));
  746. }
  747.  
  748. Poly Poly::exp(int k) const {
  749.   int lead, lz = 0;
  750.   while (lz < deg() && !a[lz])
  751.     ++lz;
  752.   if (lz == deg() && !a[lz])
  753.     return *this;
  754.   lead = a[lz];
  755.   if (lz * (ll)k > deg())
  756.     return zeroes(deg());
  757.   Poly part = Poly(vector<int>(a.begin() + lz, a.begin() + deg() - lz * (k - 1) + 1)) * nt.inv(lead);
  758.   part = (part.ln() * k).exp() * mpow(lead, k);
  759.   vector<int> ret(deg() + 1);
  760.   memcpy(ret.begin().base() + lz * k, part.base(), sizeof(int) * (deg() - lz * k + 1));
  761.   return ret;
  762. }
  763.  
  764. Poly Poly::operator/(const Poly& rhs) const {
  765.   int n = deg(), m = rhs.deg();
  766.   if (n < m)
  767.     return 0;
  768.   Poly ta(vector<int>(a.rbegin(), a.rend())),
  769.           tb(vector<int>(rhs.a.rbegin(), rhs.a.rend()));
  770.   ta.redeg(n - m);
  771.   tb.redeg(n - m);
  772.   Poly q = ta.quo(tb);
  773.   reverse(q.a.begin(), q.a.end());
  774.   return q;
  775. }
  776.  
  777. pair<Poly, Poly> Poly::div(const Poly &rhs) const {
  778.   if (deg() < rhs.deg())
  779.     return make_pair(0, *this);
  780.   int n = deg(), m = rhs.deg();
  781.   Poly q = operator/(rhs), r;
  782.   int lgn = 0;
  783.   while ((1 << lgn) < rhs.deg())
  784.     ++lgn;
  785.   int t = (1 << lgn) - 1;
  786.   r = zeroes(t);
  787.   Poly tmp = zeroes(t);
  788.   for (int i = 0; i <= m; ++i)
  789.     if ((r[i & t] += rhs[i]) >= P)
  790.       r[i & t] -= P;
  791.   for (int i = 0; i <= n - m; ++i)
  792.     if ((tmp[i & t] += q[i]) >= P)
  793.       tmp[i & t] -= P;
  794.   ntt.fft(r.base(), lgn, 1);
  795.   ntt.fft(tmp.base(), lgn, 1);
  796.   for (int i = 0; i <= t; ++i)
  797.     tmp[i] = tmp[i] * (ll)r[i] % P;
  798.   ntt.fft(tmp.base(), lgn, -1);
  799.   memset(r.base(), 0, sizeof(int) << lgn);
  800.   for (int i = 0; i <= n; ++i)
  801.     if ((r[i & t] += a[i]) >= P)
  802.       r[i & t] -= P;
  803.   for (int i = 0; i < m; ++i)
  804.     if ((r[i] -= tmp[i]) < 0)
  805.       r[i] += P;
  806.   return make_pair(q, r.slice(m - 1));
  807. }
  808.  
  809. Poly Poly::operator%(const Poly &rhs) const {
  810.   if (deg() < rhs.deg())
  811.     return *this;
  812.   return div(rhs).second;
  813. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!