首页 > 其他分享 >【模板】多项式全家桶（多项式初等函数（部分））

【模板】多项式全家桶（多项式初等函数（部分））

时间：2024-02-16 16:44:06浏览次数：35  
标签：const int 多项式 poly operator return 初等 模板 size

【模板】多项式初等函数

同时作为 https://github.com/caijianhong/template-poly 的 document。

杂项

数域为 \(\mathbb F_{998244353}\)，所以定义了 mint 为 modint<998244353>。

poly 是多项式的类型，从 std::vector<mint> 继承而来。poly 的构造函数如下：

poly();
explicit poly(int n); // n 为项数，n 个 0
poly(const vector<mint>& vec); // vec[0] 为常数项，vec[1] 为一次项，以此类推
poly(initializer_list<mint> il);

以及

  • poly& poly::cut(int lim); 效果等同于截断到 \(lim\) 项或补零至 \(lim\) 项后返回自己。

  • istream& operator>>(istream& is, poly& a); 和 ostream& operator<<(ostream& os, const poly& a); 如其名。

  • poly operator<<(poly a, const int& k); 和poly operator>>(poly a, const int& k); 分别是乘 \(x^k\) 和除 \(x^k\)（\(<k\) 次项舍弃）。有对应的 operator<<= 和 operator>>=。

  • void poly::ntt(int op); 是 NTT，\(op=1\) 是 DFT，\(op=-1\) 是 IDFT。

  • poly concalc(int n, vector<poly> vec, const function<mint(vector<mint>)>& func); 这个接口主要用于实现牛顿迭代，\(n\) 是最高次数，\(vec\) 是若干多项式，\(func\) 是一个计算的回调函数，如计算多项式乘法是这样的：

    • concalc(len, {a, b}, [](vector<mint> vec) { return vec[0] * vec[1]; });
    • 即计算 \(a\cdot b\)。\(a, b\) 都是多项式。

多项式单点求值

问题

给出有限项的多项式 \(F(x)\) 和 \(x_0\)，求 \(F(x_0)\)。

mint poly::operator()(const mint& x) const;

solution

秦九韶算法。\(O(n)\)。

多项式加法、减法、数乘

问题

给出有限项的多项式 \(F(x),G(x)\) 和 \(\lambda\)，求 \(H(x)=F(x)\pm G(x)\) 或 \(H(x)=\lambda F(x)\)。

poly operator+(poly a, const poly& b);
poly operator-(poly a, const poly& b);
poly operator*(poly a, const mint& k);
poly operator*(const mint& k, poly a);
poly operator/(poly a, const mint& k);

有对应的 operator+=，operator-=，operator*= 和 operator/=。

solution

对应位相加、相减、数乘。\(O(n)\)。

多项式乘法

问题

给出多项式 \(F(x), G(x)\)，求 \(H(x)=F(x)G(x)\)。

poly operator*(const poly& a, const poly& b);

有对应的 operator*=

solution

https://www.cnblogs.com/caijianhong/p/template-fft.html。\(O(n\log n)\)。

多项式乘法逆

问题

给出 \(F(x)\)，求 \(H(x)\bmod x^{lim}\) 满足 \(H(x)F(x)\equiv 1\pmod{x^{lim}}\)。注意，此处 \(H(x)\) 是无限项的多项式，我们只需要 \(H(x)\bmod x^{lim}\)。

poly getInv(const poly& a, int lim);

Newton's Method

给出 \(F(x), G(H(x))\)，我们需要找到 \(H(x)\) 使得 \(G(H(x))=0\)。

设 \(n\) 为偶数，已经知道了 \(H_*(x)=H(x)\bmod{x^{n/2}}\) 满足 \(G(H_*(x))\equiv 0\pmod{x^{n/2}}\)（\(H(0)\) 需要特殊计算）。想知道 \(H(x)\bmod x^{n}\)。

在 \(H(x)=H_*(x)\) 处对 \(G(H(x))\) 作泰勒展开。

\[G(H(x))=\sum_{i=0}^{+\infty}\dfrac{G^{(i)}(H_*(x))}{i!}(H(x)-H_*(x))^i=0 \]

上式，若两边 \(\bmod x^{n}\)，因为 \(H(x)-H_*(x)\) 的前 \(n/2\) 项系数全零，所以 \((H(x)-H_*(x))^i\) 在 \(i\geq 2\) 时是零。

\[G(H(x))\equiv \sum_{i=0}^{+\infty}\dfrac{G^{(i)}(H_*(x))}{i!}(H(x)-H_*(x))^i\equiv 0\pmod{x^n} \]

\[G(H(x))\equiv G(H_*(x))+G'(H_*(x))(H(x)-H_*(x))\equiv 0\pmod{x^n} \]

所以

\[G(H_*(x))+G'(H_*(x))H(x)\equiv G'(H_*(x))H_*(x)\pmod{x^n} \]

\[H(x)\equiv H_*(x)-\dfrac{G(H_*(x))}{G'(H_*(x))}\pmod {x^n} \]

注意这个 \(G'\) 是一个导数，我们最好指明它导的是 \(H(x)\) 而不是 \(x\)。这意味着与 \(H(x)\) 无关的项应视作常数。

\[H(x)\equiv H_*(x)-\dfrac{G(H_*(x))}{\dfrac{\mathrm d}{\mathrm dH(x)}G(H_*(x))}\pmod {x^n} \]

solution

需要找到 \(H(x)\)，使得 \(G(H(x))=\dfrac{1}{H(x)}-F(x)=0\)。

\[H(x)\equiv H_*(x)-\dfrac{\dfrac{1}{H_*(x)}-F(x)}{-\dfrac{1}{H_*^2(x)}}\pmod {x^n} \]

这里 \(F(x)\) 与 \(H(x)\) 无关，是常数，求导时消失了。

\[H(x)\equiv 2H_*(x)-H_*^2(x)F(x)\pmod {x^n} \]

时间复杂度 \(T(n)=T(n/2)+O(n\log n)=O(n\log n)\)。

多项式除法与取模（整除）

问题

给定一个 \(n\) 次多项式 \(F(x)\) 和一个 \(m\) 次多项式 \(G(x)\) ，请求出多项式 \(Q(x)\), \(R(x)\)，满足以下条件：

  • \(Q(x)\) 次数为 \(n-m\)，\(R(x)\) 次数小于 \(m\)（钦定，\(R(x)\) 的次数恰好为 \(m-1\)）
  • \(F(x) = Q(x) * G(x) + R(x)\)
poly operator/(poly a, poly b);
poly operator%(const poly& a, const poly& b);

有对应的 operator/= 和 operator%=。

solution

若 \(R(x)=0\)，则直接使用多项式求逆解决。不妨，定义

\[F^R(x)=x^nF\left(\dfrac{1}{x}\right) \]

注意，\(n\) 是 \(F(x)\) 最高次项的次数。发现这样 \(x^0\) 项系数成为 \(x^n\) 项系数，\(x^n\) 项系数成为 \(x^0\) 项系数，系数的顺序翻转了。而且有 \((F^R)^R(x)=F(x)\)。

\[F(1/x) = Q(1/x) * G(1/x) + R(1/x) \]

\[x^nF(1/x) = x^{n-m}Q(1/x) * x^{m}G(1/x) + x^{n-m+1}x^{m-1}R(1/x) \]

\[F^R(x)=Q^R(x)G^R(x)+x^{n-m+1}R^R(x) \]

\[F^R(x)\equiv Q^R(x)G^R(x)+x^{n-m+1}R^R(x)\pmod{x^{n-m+1}} \]

\[F^R(x)\equiv Q^R(x)G^R(x)\pmod{x^{n-m+1}} \]

消除了 \(R(x)\) 的影响，一次求逆和一次乘法可以求出 \(Q(x)\)。然后有 \(R(x)=F(x)-Q(x)G(x)\)。

\(O(n\log n)\)。

多项式形式导数、形式不定积分

问题

给出 \(F(x)\)，求 \(H(x)\) 满足 \(H(x)=\dfrac{\mathrm{d}}{\mathrm dx}F(x)\) 或者 \(H(x)=\displaystyle \int F(x)\mathrm dx\)。这里默认形式不定积分的常数项为 \(0\)，尽管应该是任意常数。

poly getDev(poly a); // 形式导数
poly getInt(poly a); // 形式不定积分

solution

\[\dfrac{\mathrm{d}}{\mathrm dx}x^k=kx^{k-1}\implies \int x^k\mathrm dx=\dfrac{1}{k+1}x^{k+1} \]

求导和不定积分都有线性性，每项分开计算。\(O(n)\)。

多项式对数函数（ln）

问题

给出 \(F(x)\)，求 \(H(x)\bmod x^{lim}\) 满足 \(H(x)\equiv \ln F(x)\pmod{x^{lim}}\)。注意，此处 \(H(x)\) 是无限项的多项式，我们只需要 \(H(x)\bmod x^{lim}\)。另外需要保证 \(F(0)=1\) 否则在这个域上不存在。

poly getLn(const poly& a, int lim);

给出 \(\ln\) 的麦克劳林级数：

\[\ln(1+x)=\sum_{i=0}^{+\infty}\dfrac{(-1)^{i+1}}{i}x^i \]

\[\ln(1-x)=-\sum_{i=0}^{+\infty}\frac{x^i}{i} \]

要点：\(\ln(1+x)\) 的关于 \(x\) 的导数是复合函数求导，是 \(\frac{1}{1+x}\)，二阶导也是复合函数，是 \(\frac{-1}{(1+x)^2}\)。

solution

复合函数求导的链式法则：\((F(G(x)))'=F'(G(x))G'(x)\) 具体来说是

\[\dfrac{\mathrm d}{\mathrm dx}F(G(x))=\dfrac{\mathrm d}{\mathrm dG(x)}F(G(x))\cdot \dfrac{\mathrm d}{\mathrm dx}G(x) \]

对 \(\ln F(x)\) 求导再积分得到

\[\dfrac{\mathrm d}{\mathrm dx}\ln F(x)=\dfrac{\dfrac{\mathrm d}{\mathrm dx}F(x)}{F(x)}\implies \ln F(x)=\int\mathrm d\ln F(x)=\int \dfrac{\dfrac{\mathrm d}{\mathrm dx}F(x)}{F(x)}\mathrm dx \]

即 \(\ln F(x)=\int F'(x)/F(x)\mathrm dx\)。注意不要把 dx 约掉。

多项式对数函数（exp）

问题

给出 \(F(x)\)，求 \(H(x)\bmod x^{lim}\) 满足 \(H(x)\equiv \exp F(x)\pmod{x^{lim}}\)。注意，此处 \(H(x)\) 是无限项的多项式，我们只需要 \(H(x)\bmod x^{lim}\)。另外需要保证 \(F(0)=0\) 否则在这个域上不存在。

poly getExp(const poly& a, int lim);

给出 \(\exp\) 的麦克劳林级数：

\[\exp x=\sum_{i=0}^{+\infty}\frac{x^i}{i!} \]

solution

需要找到 \(H(x)\)，使得 \(G(H(x))=\ln H(x)-F(x)=0\)。应用 Newton's Method。

\[H(x)\equiv H_*(x)-\dfrac{\ln H_*(x)-F(x)}{\dfrac{1}{H_*(x)}}\pmod {x^n} \]

\[H(x)\equiv H_*(x)(1-\ln H_*(x)+F(x))\pmod {x^n} \]

时间复杂度 \(T(n)=T(n/2)+O(n\log n)=O(n\log n)\)。

多项式快速幂

问题

给出 \(F(x)\) 和 \(k\)，求 \(H(x)\bmod x^{lim}\) 满足 \(H(x)\equiv F^k(x)\pmod{x^{lim}}\)。注意，此处 \(H(x)\) 是有限项但真的很多项的多项式，我们只需要 \(H(x)\bmod x^{lim}\)。

poly qpow(const poly& a, string k, int lim); // k 是高精度数
poly qpow(const poly& a, LL k, int lim);

solution

\[H(x)=F^k(x)=\exp(\ln F^k(x))=\exp(k\ln F(x)) \]

需要通过微调使得 \(F(0)=1\)，首先将 \(F(x)\) 除 \(x\) 直到有常数项，然后所有系数除掉 \(F(0)\)。最后再搞回去。这里注意的是取模问题，有两个定理：

定理一

\(p\) 是质数。来源：这里第一篇题解、这个链接暂时无法访问

\[F^p(x)\equiv F(x^p)\pmod p \]

证明略。反正可以得出，因为 \(n<p\)，所以 \(F^p(x)\equiv F(x^p)\equiv F(0)\equiv 1\pmod {x^{n}}\)（在 \(\mathbb F_p\) 上）。

所以在对 \(\ln F(x)\) 点乘 \(k\) 时，\(k\) 可以对 \(p\) 取模，这里的 \(p=998244353\)。

定理二

\(p\) 是质数，\(a\not\equiv0\pmod p\)，费马小定理：

\[a^{p-1}\equiv 1\pmod p \]

所以将常数项乘回去时，\(k\) 要对 \(p-1\) 取模，即最终的多项式要乘 \(F^{k\bmod (p-1)}(0)\) 再乘上 \(x\) 的若干次方。

多项式开根

问题

给出 \(F(x)\)，求 \(H(x)\bmod x^{lim}\) 满足 \(H^2(x)\equiv F(x)\pmod{x^{lim}}\)。注意，此处 \(H(x)\) 是无限项的多项式，我们只需要 \(H(x)\bmod x^{lim}\)。另外需要保证 \(F(0)\) 在这个域上有二次剩余。

mint sqrt(const mint& c); // 一个数求二次剩余，无解抛出异常
poly getSqrt(const poly& a, int lim);

solution

需要找到 \(H(x)\)，使得 \(G(H(x))= H^2(x)-F(x)=0\)。应用 Newton's Method。

\[H(x)\equiv H_*(x)-\dfrac{ H_*^2(x)-F(x)}{ 2H_*(x)}\pmod {x^n} \]

\[H(x)\equiv \dfrac{ H_*(x)}{ 2}+\dfrac{F(x)}{2H_*(x)}\pmod {x^n} \]

时间复杂度 \(T(n)=T(n/2)+O(n\log n)=O(n\log n)\)。

拉格朗日插值

问题

\(n\) 个点 \((x_i,y_i)\) 可以唯一地确定一个多项式 \(y = f(x)\)。现在，给定这 \(n\) 个点，请你确定这个多项式。

poly lagrange(const vector<pair<mint, mint>>& a);

solution

定义

\[\ell_j(x)=\prod_{i\neq j}\frac{x-x_i}{x_j-x_i} \]

则

\[f(x)=\sum_i\ell_i(x)y_i \]

正确性比较显然。先暴力求出 \(\prod_i(x-x_i)\)，然后每次将他 \(O(n)\) 的除掉一个二项式。\(O(n^2)\)。

Berlekamp–Massey 算法（最短线性递推式）

问题

给出一个数列 \(P\) 从 \(0\) 开始的前 \(n\) 项。求序列 \(P\) 在 \(\bmod~998244353\) 下的最短线性递推式。注意需要保证最短线性递推式长度 \(\leq n/2\)。

poly BM(poly a);

code

(solution is WIP, replace it with the code)

poly BM(poly a) {
  poly ans, lst;
  int w = 0;
  mint delta = 0;
  for (size_t i = 0; i < a.size(); i++) {
    mint tmp = -a[i];
    for (size_t j = 0; j < ans.size(); j++) tmp += ans[j] * a[i - j - 1];
    if (tmp == 0) continue;
    if (ans.empty()) {
      w = i;
      delta = tmp;
      ans = vector<mint>(i + 1, 0);
    } else {
      auto now = ans;
      mint mul = -tmp / delta;
      if (ans.size() < lst.size() + i - w) ans.resize(lst.size() + i - w);
      ans[i - w - 1] -= mul;
      for (size_t j = 0; j < lst.size(); j++) ans[i - w + j] += lst[j] * mul;
      if (now.size() <= lst.size() + i - w) {
        w = i;
        lst = now;
        delta = tmp;
      }
    }
  }
  return ans << 1;
}

Bostan-Mori 算法（求分式第 \(n\) 项）

问题

给出多项式 \(F(x), G(x)\) 和 \(k\)，求 \([x^k]F(x)/G(x)\)。

mint divide_at(poly f, poly g, T n);

code

link: https://www.cnblogs.com/Potassium/p/15130342.html

(solution is WIP, replace it with the code)

template <class T>
mint divide_at(poly f, poly g, T n) {
  for (; n; n >>= 1) {
    poly r = g;
    for (size_t i = 1; i < r.size(); i += 2) r[i] *= -1;
    f *= r;
    g *= r;
    for (size_t i = n & 1; i < f.size(); i += 2) f[i >> 1] = f[i];
    f.resize((f.size() + 1) >> 1);
    for (size_t i = 0; i < g.size(); i += 2) g[i >> 1] = g[i];
    g.resize((g.size() + 1) >> 1);
  }
  return f.empty() ? 0 : f[0] / g[0];
}

常系数齐次线性递推

问题

求一个满足 \(k\) 阶齐次线性递推数列 \({a_i}\) 的第 \(n\) 项，即：

\[a_n=\sum\limits_{i=1}^{k}f_i \times a_{n-i} \]

给出的是 \(a_0,a_1,\cdots,a_{k-1}\) 和 \(f_1, f_2, \cdots, f_k\) 注意下标。

code

link: https://www.cnblogs.com/Potassium/p/15130342.html

(solution is WIP, replace it with the code)

template <class T>
mint linear_rec(poly a, poly f, T n) {
  // a[n] = sum_i f[i] * a[n - i]
  a.resize(f.size() - 1);
  f = poly{1} - f;
  poly g = a * f;
  g.resize(a.size());
  return divide_at(g, f, n);
}

poly.cpp

#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) void(0)
#endif
typedef long long LL;
template <unsigned umod>
struct modint {
  static constexpr int mod = umod;
  unsigned v;
  modint() : v(0) {}
  template <class T, enable_if_t<is_integral<T>::value>* = nullptr>
  modint(T x) {
    x %= mod;
    if (x < 0) x += mod;
    v = x;
  }
  modint(const string& str) {
    v = 0;
    size_t i = 0;
    if (str.front() == '-') i += 1;
    while (i < str.size()) {
      assert(isdigit(str[i]));
      v = (v * 10ull % umod + str[i] - '0') % umod;
      i += 1;
    }
    if (str.front() == '-' && v) v = umod - v;
  }
  modint operator+() const { return *this; }
  modint operator-() const { return modint() - *this; }
  friend int raw(const modint& self) { return self.v; }
  friend istream& operator>>(istream& is, modint& self) {
    string str;
    is >> str;
    self = str;
    return is;
  }
  friend ostream& operator<<(ostream& os, const modint& self) {
    return os << raw(self);
  }
  modint& operator+=(const modint& rhs) {
    v += rhs.v;
    if (v >= umod) v -= umod;
    return *this;
  }
  modint& operator-=(const modint& rhs) {
    v -= rhs.v;
    if (v >= umod) v += umod;
    return *this;
  }
  modint& operator*=(const modint& rhs) {
    v = static_cast<unsigned>(1ull * v * rhs.v % umod);
    return *this;
  }
  modint& operator/=(const modint& rhs) {
    static constexpr size_t ilim = 1 << 20;
    static modint inv[ilim + 10];
    static int sz = 0;
    assert(rhs.v);
    if (rhs.v > ilim) return *this *= qpow(rhs, mod - 2);
    if (!sz) inv[1] = sz = 1;
    while (sz < (int)rhs.v) {
      for (int i = sz + 1; i <= sz << 1; i++) inv[i] = -mod / i * inv[mod % i];
      sz <<= 1;
    }
    return *this *= inv[rhs.v];
  }
  template <class T>
  friend modint qpow(modint a, T b) {
    modint r = 1;
    for (; b; b >>= 1, a *= a)
      if (b & 1) r *= a;
    return r;
  }
  friend modint operator+(modint lhs, const modint& rhs) { return lhs += rhs; }
  friend modint operator-(modint lhs, const modint& rhs) { return lhs -= rhs; }
  friend modint operator*(modint lhs, const modint& rhs) { return lhs *= rhs; }
  friend modint operator/(modint lhs, const modint& rhs) { return lhs /= rhs; }
  friend bool operator==(const modint& lhs, const modint& rhs) {
    return lhs.v == rhs.v;
  }
  friend bool operator!=(const modint& lhs, const modint& rhs) {
    return lhs.v != rhs.v;
  }
};
typedef modint<998244353> mint;
int glim(const int& x) { return 1 << (32 - __builtin_clz(x - 1)); }
int bitctz(const int& x) { return __builtin_ctz(x); }
struct poly : vector<mint> {
  poly() {}
  explicit poly(int n) : vector<mint>(n) {}
  poly(const vector<mint>& vec) : vector<mint>(vec) {}
  poly(initializer_list<mint> il) : vector<mint>(il) {}
  mint operator()(const mint& x) const;
  poly& cut(int lim);
  void ntt(int op);
};
void print(const poly& a) {
  for (size_t i = 0; i < a.size(); i++) debug("%d, ", raw(a[i]));
  debug("\n");
}
istream& operator>>(istream& is, poly& a) {
  for (auto& x : a) is >> x;
  return is;
}
ostream& operator<<(ostream& os, const poly& a) {
  bool flag = false;
  for (auto& x : a) {
    if (flag)
      os << " ";
    else
      flag = true;
    os << x;
  }
  return os;
}
mint poly::operator()(const mint& x) const {
  const auto& a = *this;
  mint res = 0;
  for (int i = (int)a.size() - 1; i >= 0; i--) {
    res = res * x + a[i];
  }
  return res;
}
poly& poly::cut(int lim) {
  resize(lim);
  return *this;
}
void poly::ntt(int op) {
  static bool wns_flag = false;
  static vector<mint> wns;
  if (!wns_flag) {
    wns_flag = true;
    for (int j = 1; j <= 23; j++) {
      wns.push_back(qpow(mint(3), raw(mint(-1)) >> j));
    }
  }
  vector<mint>& a = *this;
  int n = a.size();
  for (int i = 1, r = 0; i < n; i++) {
    r ^= n - (1 << (bitctz(n) - bitctz(i) - 1));
    if (i < r) std::swap(a[i], a[r]);
  }
  vector<mint> w(n);
  for (int k = 1, len = 2; len <= n; k <<= 1, len <<= 1) {
    mint wn = wns[bitctz(k)];
    for (int i = raw(w[0] = 1); i < k; i++) w[i] = w[i - 1] * wn;
    for (int i = 0; i < n; i += len) {
      for (int j = 0; j < k; j++) {
        mint x = a[i + j], y = a[i + j + k] * w[j];
        a[i + j] = x + y, a[i + j + k] = x - y;
      }
    }
  }
  if (op == -1) {
    mint iz = mint(1) / n;
    for (int i = 0; i < n; i++) a[i] *= iz;
    reverse(a.begin() + 1, a.end());
  }
}
poly concalc(int n, vector<poly> vec,
             const function<mint(vector<mint>)>& func) {
  int lim = glim(n);
  int m = vec.size();
  for (auto& f : vec) f.resize(lim), f.ntt(1);
  vector<mint> tmp(m);
  poly ret(lim);
  for (int i = 0; i < lim; i++) {
    for (int j = 0; j < m; j++) tmp[j] = vec[j][i];
    ret[i] = func(tmp);
  }
  ret.ntt(-1);
  return ret;
}
poly getInv(const poly& a, int lim) {
  poly b{1 / a[0]};
  for (int len = 2; len <= glim(lim); len <<= 1) {
    poly c = vector<mint>(a.begin(), a.begin() + min(len, (int)a.size()));
    b = concalc(len << 1, {b, c}, [](vector<mint> vec) {
          return vec[0] * (2 - vec[0] * vec[1]);
        }).cut(len);
  }
  return b.cut(lim);
}
poly operator+=(poly& a, const poly& b) {
  if (a.size() < b.size()) a.resize(b.size());
  for (size_t i = 0; i < b.size(); i++) a[i] += b[i];
  return a;
}
poly operator-=(poly& a, const poly& b) {
  if (a.size() < b.size()) a.resize(b.size());
  for (size_t i = 0; i < b.size(); i++) a[i] -= b[i];
  return a;
}
poly operator*=(poly& a, const mint& k) {
  if (k == 1) return a;
  for (size_t i = 0; i < a.size(); i++) a[i] *= k;
  return a;
}
poly operator/=(poly& a, const mint& k) { return a *= 1 / k; }
poly operator<<=(poly& a, const int& k) {
  // mnltiple by x^k
  a.insert(a.begin(), k, 0);
  return a;
}
poly operator>>=(poly& a, const int& k) {
  // divide by x^k
  a.erase(a.begin(), a.begin() + min(k, (int)a.size()));
  return a;
}
poly operator*(const poly& a, const poly& b) {
  if (a.empty() || b.empty()) return {};
  int rlen = a.size() + b.size() - 1;
  int len = glim(rlen);
  if (1ull * a.size() * b.size() <= 1ull * len * bitctz(len)) {
    poly ret(rlen);
    for (size_t i = 0; i < a.size(); i++)
      for (size_t j = 0; j < b.size(); j++) ret[i + j] += a[i] * b[j];
    return ret;
  } else {
    return concalc(len, {a, b},
                   [](vector<mint> vec) { return vec[0] * vec[1]; })
        .cut(rlen);
  }
}
poly operator/(poly a, poly b) {
  if (a.size() < b.size()) return {};
  int rlen = a.size() - b.size() + 1;
  reverse(a.begin(), a.end());
  reverse(b.begin(), b.end());
  a = (a * getInv(b, rlen)).cut(rlen);
  reverse(a.begin(), a.end());
  return a;
}
poly operator-(poly a, const poly& b) { return a -= b; }
poly operator%(const poly& a, const poly& b) {
  return (a - (a / b) * b).cut(b.size() - 1);
}
poly operator*=(poly& a, const poly& b) { return a = a * b; }
poly operator/=(poly& a, const poly& b) { return a = a / b; }
poly operator%=(poly& a, const poly& b) { return a = a % b; }
poly operator+(poly a, const poly& b) { return a += b; }
poly operator*(poly a, const mint& k) { return a *= k; }
poly operator*(const mint& k, poly a) { return a *= k; }
poly operator/(poly a, const mint& k) { return a /= k; }
poly operator<<(poly a, const int& k) { return a <<= k; }
poly operator>>(poly a, const int& k) { return a >>= k; }
poly getDev(poly a) {
  a >>= 1;
  for (size_t i = 1; i < a.size(); i++) a[i] *= i + 1;
  return a;
}
poly getInt(poly a) {
  a <<= 1;
  for (size_t i = 1; i < a.size(); i++) a[i] /= i;
  return a;
}
poly getLn(const poly& a, int lim) {
  assert(a[0] == 1);
  return getInt(getDev(a) * getInv(a, lim)).cut(lim);
}
poly getExp(const poly& a, int lim) {
  assert(a[0] == 0);
  poly b{1};
  for (int len = 2; len <= glim(lim); len <<= 1) {
    poly c = vector<mint>(a.begin(), a.begin() + min(len, (int)a.size()));
    b = concalc(len << 1, {b, getLn(b, len), c}, [](vector<mint> vec) {
          return vec[0] * (1 - vec[1] + vec[2]);
        }).cut(len);
  }
  return b.cut(lim);
}
poly qpow(const poly& a, string k, int lim) {
  size_t i = 0;
  while (i < a.size() && a[i] == 0) i += 1;
  if (i == a.size() || (i > 0 && k.size() >= 9) ||
      1ull * i * raw(mint(k)) >= 1ull * lim)
    return poly(lim);
  lim -= i * raw(mint(k));
  return getExp(getLn(a / a[i] >> i, lim) * k, lim) *
             qpow(a[i], raw(modint<mint::mod - 1>(k)))
         << i * raw(mint(k));
}
mint sqrt(const mint& c) {
  static const auto check = [](mint c) {
    return qpow(c, (mint::mod - 1) >> 1) == 1;
  };
  if (raw(c) <= 1) return 1;
  if (!check(c)) throw "No solution!";
  static mt19937 rng{random_device{}()};
  mint a = rng();
  while (check(a * a - c)) a = rng();
  typedef pair<mint, mint> number;
  const auto mul = [=](number x, number y) {
    return make_pair(x.first * y.first + x.second * y.second * (a * a - c),
                     x.first * y.second + x.second * y.first);
  };
  const auto qpow = [=](number a, int b) {
    number r = {1, 0};
    for (; b; b >>= 1, a = mul(a, a))
      if (b & 1) r = mul(r, a);
    return r;
  };
  mint ret = qpow({a, 1}, (mint::mod + 1) >> 1).first;
  return min(raw(ret), raw(-ret));
}
poly getSqrt(const poly& a, int lim) {
  poly b{sqrt(a[0])};
  for (int len = 2; len <= glim(lim); len <<= 1) {
    poly c = vector<mint>(a.begin(), a.begin() + min(len, (int)a.size()));
    b = (c * getInv(b * 2, len) + b / 2).cut(len);
  }
  return b.cut(lim);
}
template <class T>
mint divide_at(poly f, poly g, T n) {
  for (; n; n >>= 1) {
    poly r = g;
    for (size_t i = 1; i < r.size(); i += 2) r[i] *= -1;
    f *= r;
    g *= r;
    for (size_t i = n & 1; i < f.size(); i += 2) f[i >> 1] = f[i];
    f.resize((f.size() + 1) >> 1);
    for (size_t i = 0; i < g.size(); i += 2) g[i >> 1] = g[i];
    g.resize((g.size() + 1) >> 1);
  }
  return f.empty() ? 0 : f[0] / g[0];
}
template <class T>
mint linear_rec(poly a, poly f, T n) {
  // a[n] = sum_i f[i] * a[n - i]
  a.resize(f.size() - 1);
  f = poly{1} - f;
  poly g = a * f;
  g.resize(a.size());
  return divide_at(g, f, n);
}
poly BM(poly a) {
  poly ans, lst;
  int w = 0;
  mint delta = 0;
  for (size_t i = 0; i < a.size(); i++) {
    mint tmp = -a[i];
    for (size_t j = 0; j < ans.size(); j++) tmp += ans[j] * a[i - j - 1];
    if (tmp == 0) continue;
    if (ans.empty()) {
      w = i;
      delta = tmp;
      ans = vector<mint>(i + 1, 0);
    } else {
      auto now = ans;
      mint mul = -tmp / delta;
      if (ans.size() < lst.size() + i - w) ans.resize(lst.size() + i - w);
      ans[i - w - 1] -= mul;
      for (size_t j = 0; j < lst.size(); j++) ans[i - w + j] += lst[j] * mul;
      if (now.size() <= lst.size() + i - w) {
        w = i;
        lst = now;
        delta = tmp;
      }
    }
  }
  return ans << 1;
}
poly lagrange(const vector<pair<mint, mint>>& a) {
  poly ans(a.size()), product{1};
  for (size_t i = 0; i < a.size(); i++) {
    product *= poly{-a[i].first, 1};
  }
  auto divide2 = [&](poly a, mint b) {
    poly res(a.size() - 1);
    for (size_t i = (int)a.size() - 1; i >= 1; i--) {
      res[i - 1] = a[i];
      a[i - 1] -= a[i] * b;
    }
    return res;
  };
  for (size_t i = 0; i < a.size(); i++) {
    mint denos = 1;
    for (size_t j = 0; j < a.size(); j++) {
      if (i != j) denos *= a[i].first - a[j].first;
    }
    poly numes = divide2(product, -a[i].first);
    ans += a[i].second / denos * numes;
  }
  return ans;
}
int main() {
#ifndef LOCAL
  cin.tie(nullptr)->sync_with_stdio(false);
#endif
  int n;
  string k;
  cin >> n >> k;
  poly a(n + 1);
  cin >> a;
  cout << getDev(qpow(getLn(a + poly{2 - a[0]} -
                                getExp(getInt(getInv(getSqrt(a, n), n)), n + 1),
                            n + 1) +
                          poly{1},
                      k, n + 1))
       << endl;
}

标签：const,int,多项式,poly,operator,return,初等,模板,size
From： https://www.cnblogs.com/caijianhong/p/18017272

相关文章

  • 【题解】P4722 【模板】最大流 加强版 / 预流推进
    更好阅读体验CHAPTER0废话1.常见的最大流算法可以大致分为两种：找增广路和预流推进2.从理论时间复杂度上界来看，找增广路算法的最坏时间复杂度为\(O(n^2m)\)，而预流推进算法的最坏时间复杂度为\(O(n^2\sqrt{m})\)，看起来预流推进要显然优于找增广路，但在实际应用（尤指OI）中，由于包......
  • 多项式从入门到进门
    多项式全家福多项式一个以\(x\)为变量的多项式定义在一个代数域\(F\)上，将函数\(A(x)\)表示为形式和：\[A(x)=\sum\limits_{i=0}^{n-1}a_ix^i\]我们称\(a_0,a_1,\cdots,a_n\)为多项式的系数，所有系数都属于数域\(\mathbbF\)。如果一个多项式的最高次的非零系数是\(k......
  • 归并排序模板
    #include<bits/stdc++.h>usingnamespacestd;constintN=1e5+10;intn,s[N],res[N];voidmerge_sort(ints[],intl,intr){intmid=(l+r)>>1;if(l>=r)return;merge_sort(s,l,mid);merge_sort(s,mid+1,r);inti=l,k=0,j......
  • 二叉树遍历问题模板
    在二叉树遍历问题中，有三种常见的遍历方式：前序遍历、中序遍历和后序遍历。以下是这三种遍历方式的递归模板：1.前序遍历（PreorderTraversal）：defpreorderTraversal(root):ifnotroot:return[]result=[]result.append(root.val)#处理当前节点......
  • 回溯算法模板
    回溯算法的模板通常包含递归函数和回溯过程。以下是一个通用的回溯算法模板：defbacktrack(start,path,other_parameters):#满足结束条件时，将当前路径加入结果ifsatisfies_end_condition:result.append(path[:])return#从start开始遍历可......
  • 模板
    01.模版的概念（了解）1.函数或类是通用，但是里面的数据类型的多种状态2.模版有：函数和类02.函数模版（重点）1.什么是函数模版函数模板，实际上是建立一个通用函数，其函数类型和形参类型不具体制定，用一个虚拟的类型来代表。这个通用函数就成为函数模板2.怎么编写函数模版//T代表泛型的......
  • P3367 【模板】并查集
    原题链接并查集模板练手。递归版本#include<bits/stdc++.h>usingnamespacestd;constintN=1e4+5;intfather[N];intfind(intmid){if(father[mid]!=mid){father[mid]=find(father[mid]);}returnfather[mid];}voidunion_fa(intx,inty){......
  • P4512 【模板】多项式除法
    为什么不能直接用\(F(x)\timesG(x)^{-1}\)做？\(G(x)^{-1}\)是模\(x^{m+1}\)意义下的，\(n-m>m\)时得到的\(Q(x)\)就是错的。记\(F'(x)\)为\(F(x)\)系数翻转后的多项式，即若\(F(x)=\sum\limits_{i=0}^nf_ix^i\)，则\(F'(x)=\sum\limits_{i=0}^nf_{n......
  • P3811 【模板】模意义下的乘法逆元
    原题链接题解由于时间限制过于严苛，遂采用线性递推方式\(p=k·i+b\)，\((1\leqslantb<r<p)\)\(k·i+b=0\)\((mod\p)\)同时乘上\(i^{-1}\b^{-1}\)\(k·b^{-1}+i^{-1}=0\(mod\p)\)\(i^{-1}=-k·b^{-1}\(mod\p)\)\(i^{-1}=(-[\frac{p}{i}]+p)+(p\mod\i)^{-1}......
  • 各大排序的模板
    1.冒泡排序 1for(i=n;i>=1;--i)2{3for(j=1;j<=i;++j)4{5if(a[j]>a[j+1])6{7swap(a[j],a[j+1]);8}9}10}   2.快速排序1.懒人函数 1sort(a+1,a+n+1);   2.正常的1vo......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99