跳转至

多项式初等函数

本文转载(或修改)自 OI-Wiki

本页面包含多项式常见的初等函数操作。具体而言,本页面包含如下内容:

  1. 多项式求逆
  2. 多项式开方
  3. 多项式除法
  4. 多项式取模
  5. 多项式指数函数
  6. 多项式对数函数
  7. 多项式快速幂
  8. 多项式三角函数
  9. 多项式反三角函数
初等函数与非初等函数

初等函数的定义如下:

若域 \(F\) 中存在映射 \(u\to \partial u\) 满足:

  1. \(\partial(u+v)=\partial u+\partial v\)
  2. \(\partial(uv)=u\partial v+v\partial u\)

则称这个域为 微分域

若微分域 \(F\) 上的函数 \(u\) 满足以下的任意一条条件,则称该函数 \(u\) 为初等函数:

  1. \(u\)\(F\) 上的代数函数。
  2. \(u\)\(F\) 上的指数性函数,即存在 \(a\in F\) 使得 \(\partial u=u\partial a\) .
  3. \(u\)\(F\) 上的对数性函数,即存在 \(a\in F\) 使得 \(\partial u=\frac{\partial a}{a}\) .

以下是常见的初等函数:

  1. 代数函数:存在有限次多项式 \(P\) 使得 \(P(f(x))=0\) 的函数 \(f(x)\),如 \(2x+1\) , \(\sqrt{x}\) , \((1+x^2)^{-1}\) , \(|x|\) .
  2. 指数函数
  3. 对数函数
  4. 三角函数
  5. 反三角函数
  6. 双曲函数
  7. 反双曲函数
  8. 以上函数的复合,如:

    \[ \frac{\mathrm{e}^{\tan x}}{1+x^2}\sin\left(\sqrt{1+\ln^2 x}\right) \]
    \[ -\mathrm{i} \ln\left(x+\mathrm{i}\sqrt{1-x^2}\right) \]

以下是常见的非初等函数:

  1. 误差函数:

    \[ \operatorname{erf}(x):=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp\left(-t^2\right)\mathrm{d}t \]

多项式求逆

P4238 【模板】多项式乘法逆

给定一个多项式 \(F(x)\) ,请求出一个多项式 \(G(x)\), 满足 \(F(x) * G(x) \equiv 1 \pmod{x^n}\)

系数对 \(998244353\) 取模。

对于 \(100\%\) 的数据,\(1 \leq n \leq 10^5\)\(0 \leq a_i \leq 10^9\)

首先,易知

\[ \left[x^{0}\right]f^{-1}\left(x\right)=\left(\left[x^{0}\right]f\left(x\right)\right)^{-1} \]

假设现在已经求出了 \(f\left(x\right)\) 在模 \(x^{\left\lceil\frac{n}{2}\right\rceil}\) 意义下的逆元 \(f^{-1}_{0}\left(x\right)\)。 有:

\[ \begin{aligned} f\left(x\right)f^{-1}_{0}\left(x\right)&\equiv 1 &\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\\ f\left(x\right)f^{-1}\left(x\right)&\equiv 1 &\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\\ f^{-1}\left(x\right)-f^{-1}_{0}\left(x\right)&\equiv 0 &\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} \end{aligned} \]

两边平方可得:

\[ f^{-2}\left(x\right)-2f^{-1}\left(x\right)f^{-1}_{0}\left(x\right)+f^{-2}_{0}\left(x\right)\equiv 0 \pmod{x^{n}} \]

两边同乘 \(f\left(x\right)\) 并移项可得:

\[ f^{-1}\left(x\right)\equiv f^{-1}_{0}\left(x\right)\left(2-f\left(x\right)f^{-1}_{0}\left(x\right)\right) \pmod{x^{n}} \]

递归计算即可,初始边界是 \(\left[x^{0}\right]f^{-1}\left(x\right)=\left(\left[x^{0}\right]f\left(x\right)\right)^{-1}\),不过改成自底而上的迭代会好得多。

时间复杂度

\[ T\left(n\right)=T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log{n}\right) \]
实现(迭代版)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
void Inv(int n, ll *F, ll *G)
{
    static ll A[MAXN], B[MAXN];
    G[0] = inv(F[0]);
    int len, lim;
    for (len = 1; len < (n << 1); len <<= 1)
    {
        lim = len << 1;
        copy(F, F + len, A), copy(G, G + len, B);
        for (int i = 0; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * len);
        NTT(A, lim, 1), NTT(B, lim, 1);
        for (int i = 0; i < lim; i++)
            G[i] = ((2 - A[i] * B[i] % p) * B[i] % p + p) % p;
        NTT(G, lim, -1);
        fill(G + len, G + lim, 0);
    }
    fill(A, A + len, 0), fill(B, B + len, 0), fill(G + n, G + len, 0);
}

多项式开方

P5205 【模板】多项式开根 P5277 【模板】多项式开根(加强版)

给定多项式 \(f\left(x\right)\),求 \(g\left(x\right)\),满足:

\[ g^{2}\left(x\right)\equiv f\left(x\right) \pmod{x^{n}} \]

普通版保证 \(f_0 = 1\),加强版只保证 \(f_0\) 有二次剩余。

首先讨论 \(\left[x^0\right]f(x)\) 不为 \(0\) 的情况。

易知:

\[ \left[x^0\right]g(x) = \sqrt{\left[x^0\right]f(x)} \]

\(\left[x^0\right]f(x)\) 没有平方根,则多项式 \(f(x)\) 没有平方根。

\(\left[x^0\right]f(x)\) 可能有多个平方根,选取不同的根会求出不同的 \(g(x)\)

假设现在已经求出了 \(f\left(x\right)\) 在模 \(x^{\left\lceil\frac{n}{2}\right\rceil}\) 意义下的平方根 \(g_{0}\left(x\right)\),则有:

\[ \begin{aligned} g_{0}^{2}\left(x\right)&\equiv f\left(x\right) &\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\\ g_{0}^{2}\left(x\right)-f\left(x\right)&\equiv 0 &\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\\ \left(g_{0}^{2}\left(x\right)-f\left(x\right)\right)^{2}&\equiv 0 &\pmod{x^{n}}\\ \left(g_{0}^{2}\left(x\right)+f\left(x\right)\right)^{2}&\equiv 4g_{0}^{2}\left(x\right)f\left(x\right) &\pmod{x^{n}}\\ \left(\frac{g_{0}^{2}\left(x\right)+f\left(x\right)}{2g_{0}\left(x\right)}\right)^{2}&\equiv f\left(x\right) &\pmod{x^{n}}\\ \frac{g_{0}^{2}\left(x\right)+f\left(x\right)}{2g_{0}\left(x\right)}&\equiv g\left(x\right) &\pmod{x^{n}}\\ 2^{-1}g_{0}\left(x\right)+2^{-1}g_{0}^{-1}\left(x\right)f\left(x\right)&\equiv g\left(x\right) &\pmod{x^{n}} \end{aligned} \]

倍增计算即可。

时间复杂度

\[ T\left(n\right)=T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log{n}\right) \]

还有一种常数较小的写法就是在倍增维护 \(g\left(x\right)\) 的时候同时维护 \(g^{-1}\left(x\right)\) 而不是每次都求逆。

\(\left[x^{0}\right]f\left(x\right)\neq 1\) 时,可能需要使用二次剩余来计算 \(\left[x^{0}\right]g\left(x\right)\)

上述方法需要知道 \(g_{0}(x)\) 的逆,所以常数项不能为 \(0\)

\(\left[x^0\right]f(x) = 0\),则将 \(f(x)\) 分解成 \(x^{k}h(x)\),其中 \(\left[x^0\right]h(x) \not = 0\)

  • \(k\) 是奇数,则 \(f(x)\) 没有平方根。
  • \(k\) 是偶数,则求出 \(h(x)\) 的平方根 \(\sqrt{h(x)}\),然后得到 \(g(x) \equiv x^{k/2} \sqrt{h(x)} \pmod{x^{n}}\)
实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
void Sqrt(ll *F, ll *G, int n) // G = sqrt(F)
{
    static const ll inv2 = inv(2);
    G[0] = 1;
    static ll A[MAXN], B[MAXN];
    int len, lim;
    for (len = 1; len < (n << 1); len <<= 1)
    {
        lim = len << 1;
        copy(F, F + len, A);
        Inv(G, B, len);
        for (int i = 0; i < lim; i++)
            rev[i] = (rev[i >> 1] >> 1) | ((lim >> 1) * (i & 1));
        NTT(A, lim, 1), NTT(B, lim, 1);
        for (int i = 0; i < lim; i++) A[i] = A[i] * B[i] % p;
        NTT(A, lim, -1);
        for (int i = 0; i < len; i++) G[i] = (G[i] + A[i]) % p * inv2 % p;
        fill(G + len, G + lim, 0);
    }
    fill(A, A + len, 0), fill(B, B + len, 0), fill(G + n, G + len, 0);
}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void Sqrt(ll *F, ll *G, int n)
{
    static const ll inv2 = inv(2);
    ll x1, x2;
    cipolla::cipolla(F[0], x1, x2);
    G[0] = min(x1, x2);
    static ll A[MAXN], B[MAXN];
    int len, lim;
    for (len = 1; len < (n << 1); len <<= 1)
    {
        lim = len << 1;
        copy(F, F + len, A);
        Inv(G, B, len);
        for (int i = 0; i < lim; i++)
            rev[i] = (rev[i >> 1] >> 1) | ((lim >> 1) * (i & 1));
        NTT(A, lim, 1), NTT(B, lim, 1);
        for (int i = 0; i < lim; i++) A[i] = A[i] * B[i] % p;
        NTT(A, lim, -1);
        for (int i = 0; i < len; i++) G[i] = (G[i] + A[i]) % p * inv2 % p;
        fill(G + len, G + lim, 0);
    }
    fill(A, A + len, 0), fill(B, B + len, 0), fill(G + n, G + len, 0);
}

多项式除法 & 取模

P4512 【模板】多项式除法

给定多项式 \(F\left(x\right),G\left(x\right)\),求 \(G\left(x\right) \div F\left(x\right)\) 的商式 \(Q\left(x\right)\) 和余式 \(R\left(x\right)\)

发现若能消除 \(R\left(x\right)\) 的影响则可直接 多项式求逆 解决。

具体来说,设多项式 \(A\)\(n\) 次多项式,考虑一种操作 \(R\),使得

\[A_R(x) = x^n A\left(\frac{1}{x}\right)\]

稍微想象一下,可以发现 \(A_R[i] = A[n - i]\)\([i]\) 表示多项式的第 \(i\) 次系数)。这个操作可以 \(O(n)\) 完成。

然后开始推式子:

\[\begin{aligned} F(x) &= Q(x) * G(x) + R(x) \\ F\left(\frac{1}{x}\right) &= Q\left(\frac{1}{x}\right) * G\left(\frac{1}{x}\right) + R\left(\frac{1}{x}\right) \\ x^n F\left(\frac{1}{x}\right) &= x^{n - m} Q\left(\frac{1}{x}\right) * x^m G\left(\frac{1}{x}\right) + x^{n - m + 1} * x^{m - 1} R\left(\frac{1}{x}\right) \\ 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}} \\ Q_R(x) &\equiv F_R(x) * G_R^{-1}(x) \pmod{x^{n-m+1}} \end{aligned}\]

求一遍 \(G_R\) 的逆,然后就可以利用多项式乘法求出 \(Q\)。然后

\[R(x) = F(x) - G(x) * Q(x)\]

直接计算即可。时间复杂度 \(O(n \log n)\)

时间复杂度 \(O\left(n\log{n}\right)\)

实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
void Div(ll *F, ll *G, ll *Q, ll *R, int n, int m)
{
    static ll B[MAXN], BR[MAXN], A[MAXN];
    static int lim;
    for (lim = 1; lim < (n << 1); lim <<= 1);
    copy(G, G + m, B), reverse(B, B + m), Inv(n, B, BR);
    for (lim = 1; lim < (n << 1); lim <<= 1);
    copy(F, F + n, B), reverse(B, B + n);
    NTT(B, lim, 1), NTT(BR, lim, 1);
    for (int i = 0; i < lim; i++) Q[i] = B[i] * BR[i] % p;
    NTT(Q, lim, -1), fill(Q + n - m + 1, Q + lim, 0), reverse(Q, Q + n - m + 1),
        fill(B, B + lim, 0), fill(BR, BR + lim, 0), copy(Q, Q + n - m + 1, A),
        copy(G, G + m, B), NTT(A, lim, 1), NTT(B, lim, 1);
    for (int i = 0; i < lim; i++) A[i] = A[i] * B[i] % p;
    NTT(A, lim, -1);
    for (int i = 0; i < m - 1; i++) R[i] = (F[i] + p - A[i]) % p;
    fill(B, B + lim, 0), fill(BR, BR + lim, 0), fill(A, A + lim, 0);
}

多项式对数函数

P4725 【模板】多项式对数函数(多项式 ln)

给定多项式 \(f(x)\),求模 \(x^{n}\) 意义下的 \(\ln{f(x)}\)

首先,对于多项式 \(f(x)\),若 \(\ln{f(x)}\) 存在,则由其定义,其必须满足:

\[ [x^{0}]f(x)=1 \]

\(\ln{f(x)}\) 求导再积分,可得:

\[ \begin{aligned} \frac{\mathrm{d} \ln{f(x)}}{\mathrm{d} x} & \equiv \frac{f'(x)}{f(x)} & \pmod{x^{n}} \\ \ln{f(x)} & \equiv \int\frac{f'(x)}{f(x)} \mathrm{d} x & \pmod{x^{n}} \end{aligned} \]

多项式的求导,积分时间复杂度为 \(O(n)\),求逆时间复杂度为 \(O(n\log{n})\),故多项式求 \(\ln\) 时间复杂度 \(O(n\log{n})\)

实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
void Ln(int n, ll *F, ll *G)
{
    static ll A[MAXN], B[MAXN];
    Inv(n, F, A), Dert(n, F, B);
    int lim;
    for (lim = 1; lim < (n << 1); lim <<= 1);
    for (int i = 0; i < lim; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (lim >> 1));
    NTT(A, lim, 1), NTT(B, lim, 1);
    for (int i = 0; i < lim; i++) A[i] = A[i] * B[i] % p;
    NTT(A, lim, -1), fill(A + n, A + lim, 0), Int(n - 1, A, G), fill(A, A + lim, 0), fill(B, B + lim, 0);
}

多项式指数函数

P4726 【模板】多项式指数函数(多项式 exp)

给定多项式 \(f(x)\),求模 \(x^{n}\) 意义下的 \(\exp{f(x)}\)

首先,对于多项式 \(f(x)\),若 \(\exp{f(x)}\) 存在,则其必须满足:

\[ [x^{0}]f(x)=0 \]

否则 \(\exp{f(x)}\) 的常数项不收敛。

\(\exp{f(x)}\) 求导,可得:

\[ \frac{\mathrm{d} \exp{f(x)}}{\mathrm{d} x} \equiv \exp{f(x)}f'(x)\pmod{x^{n}} \]

比较两边系数可得:

\[ [x^{n-1}]\frac{\mathrm{d} \exp{f(x)}}{\mathrm{d} x} = \sum_{i = 0}^{n - 1} \left([x^{i}]\exp{f(x)}\right) \left([x^{n-i-1}]f'(x)\right) \]
\[ n[x^{n}]\exp{f(x)} = \sum_{i = 0}^{n - 1} \left([x^{i}]\exp{f(x)}\right) \left((n - i)[x^{n - i}]f(x)\right) \]

使用分治 FFT 即可解决。

时间复杂度 \(O(n\log^{2}{n})\)

更加常见且快速的做法是应用 Newton's Method,时间复杂度 \(O(n\log n)\),下面是其实现:

实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
void Exp(int n, ll *F, ll *G)
{
    static ll A[MAXN], B[MAXN];
    int lim, len;
    G[0] = 1;
    for (len = 1; len < (n << 1); len <<= 1)
    {
        lim = len << 1;
        copy(F, F + len, A), Ln(len, G, B);
        for (int i = 0; i < lim; i++)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * len);
        NTT(A, lim, 1), NTT(B, lim, 1), NTT(G, lim, 1);
        for (int i = 0; i < lim; i++)
            G[i] = (1 - B[i] + A[i] + p) % p * G[i] % p;
        NTT(G, lim, -1), fill(A, A + lim, 0), fill(B, B + lim, 0);
    }
    fill(G + n, G + len, 0);
}

多项式快速幂

P5245 【模板】多项式快速幂 P5273 【模板】多项式快速幂(加强版)

计算 \(f^{k}(x)\)

普通做法为直接套用实数域中的快速幂,时间复杂度 \(O(n\log{n}\log{k})\),不够优秀,Luogu上的板子中,\(k = 10^{10^{15}}\)

(普通版)当 \([x^{0}]f(x)=1\) 时,有:

\[f^{k}(x)=\exp{\left(k\ln{f(x)}\right)}\]

(加强版)当 \([x^{0}]f(x)\neq 1\) 时,设 \(f(x)\) 的最低次项为 \(f_{i}x^{i}\),则:

\[f^{k}(x)=f_{i}^{k}x^{ik}\exp{\left(k\ln{\frac{f(x)}{f_{i}x^{i}}}\right)}\]

时间复杂度 \(O(n\log{n})\)

实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void Pow(int n, ll k, ll *F, ll *G)
{
    static ll A[MAXN], B[MAXN];
    copy(F, F + n, A);
    Ln(n, A, B);
    for (int i = 0; i < n; i++)
        B[i] = B[i] * k % p;
    Exp(n, B, G), fill(A, A + n, 0), fill(B, B + n, 0);
}
int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    initinv(4e5);
    cin >> n >> sk;
    for (char ch : sk)
        k = (k * 10 % p + (ch ^ 48)) % p;
    for (int i = 0; i < n; i++)
        cin >> F[i];
    Pow(n, k, F, G);
    for (int i = 0; i < n; i++)
        cout << G[i] << " ";
    return 0;
}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
void Pow(ll n, ll k, ll k2, ll *F, ll *G)
{
    static ll A[MAXN], B[MAXN];
    ll t = 0;
    while (t < n && !F[t]) t++;
    copy(F + t, F + n, A);
    ll m = n - t;
    ll f0 = A[0], invf0 = getinv(f0), pf0 = qpow(f0, k2);
    for (ll i = 0; i < m; i++) A[i] = A[i] * invf0 % p;
    Ln(m, A, B);
    for (ll i = 0; i < m; i++) B[i] = B[i] * k % p;
    fill(A, A + m, 0), Exp(m, B, A), t *= k2;
    for (ll i = n - t - 1; i >= 0; i--) G[i + t] = A[i] * pf0 % p;
    fill(A, A + m, 0), fill(B, B + m, 0);
}
int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    initinv(4e6);
    cin >> n >> sk;
    for (char ch : sk)
        k = (k * 10 % p + (ch ^ 48)) % p,
        k2 = (k2 * 10 % (p - 1) + (ch ^ 48)) % (p - 1);
    for (int i = 0; i < n; i++) cin >> F[i];
    if (!(F[0] == 0 && sk.size() > log10(n) + 1)) Pow(n, k, k2, F, G);
    for (int i = 0; i < n; i++) cout << G[i] << " ";
    return 0;
}

多项式三角函数

P5264 多项式三角函数

给定多项式 \(f\left(x\right)\),求模 \(x^{n}\) 意义下的 \(\sin{f\left(x\right)}\)\(\cos{f\left(x\right)}\)

首先由欧拉公式 \(\left(\mathrm{e}^{\mathrm{i}x} = \cos{x} + \mathrm{i}\sin{x}\right)\) 可以得到三角函数的另一个表达式

\[ \begin{aligned} \sin{x} &= \frac{\mathrm{e}^{\mathrm{i}x} - \mathrm{e}^{-\mathrm{i}x}}{2\mathrm{i}} \\ \cos{x} &= \frac{\mathrm{e}^{\mathrm{i}x} + \mathrm{e}^{-\mathrm{i}x}}{2} \end{aligned} \]

那么代入 \(f\left(x\right)\) 就有:

\[ \begin{aligned} \sin{f\left(x\right)} &= \frac{\exp{\left(\mathrm{i}f\left(x\right)\right)} - \exp{\left(-\mathrm{i}f\left(x\right)\right)}}{2\mathrm{i}} \\ \cos{f\left(x\right)} &= \frac{\exp{\left(\mathrm{i}f\left(x\right)\right)} + \exp{\left(-\mathrm{i}f\left(x\right)\right)}}{2} \end{aligned} \]

直接按上述表达式编写程序即可得到模 \(x^{n}\) 意义下的 \(\sin{f\left(x\right)}\)\(\cos{f\left(x\right)}\)。再由 \(\tan{f\left(x\right)} = \frac{\sin{f\left(x\right)}}{\cos{f\left(x\right)}}\) 可求得 \(\tan{f\left(x\right)}\)

模意义下的虚数

注意到我们是在 \(\mathbb{Z}_{998244353}\) 上做 NTT,那么相应地,虚数单位 \(\mathrm{i}\) 应该被换成 \(86583718\) (或 \(911660635\)):

\[ \begin{aligned} & \mathrm{i} = \sqrt{-1} \equiv \sqrt{998244352} \pmod{998244353} \\ \implies & \phantom{\text{or}} \quad \mathrm{i} \equiv 86583718 \pmod{998244353} \end{aligned} \]

对应 \(g^{\frac{p - 1}{4}}\)

实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
const ll Imag = getinv((p - 1) / 4);
void Sin(int n, ll *F, ll *G)
{
    static ll A[MAXN], B[MAXN], H[MAXN];
    for (int i = 0; i < n; i++) A[i] = F[i] * Imag;
    Exp(n, A, B), Inv(n, B, H);
    static ll inv2imag = getinv(2 * Imag);
    for (int i = 0; i < n; i++) G[i] = (B[i] - H[i] + p) % p * inv2imag % p;
    fill(A, A + n, 0), fill(B, B + n, 0), fill(H, H + n, 0);
}
void Cos(int n, ll *F, ll *G)
{
    static ll A[MAXN], B[MAXN], H[MAXN];
    for (int i = 0; i < n; i++) A[i] = F[i] * Imag;
    Exp(n, A, B), Inv(n, B, H);
    static ll inv2 = getinv(2);
    for (int i = 0; i < n; i++) G[i] = (B[i] + H[i]) % p * inv2 % p;
    fill(A, A + n, 0), fill(B, B + n, 0), fill(H, H + n, 0);
}

多项式反三角函数

P5265 多项式反三角函数

给定多项式 \(f\left(x\right)\),求模 \(x^{n}\) 意义下的 \(\arcsin{f\left(x\right)}\)\(\arctan{f\left(x\right)}\)

仿照求多项式 \(\ln\) 的方法,对反三角函数求导再积分可得:

\[ \begin{aligned} \frac{\mathrm{d}}{\mathrm{d} x} \arcsin{x} &= \frac{1}{\sqrt{1 - x^{2}}} \\ \arcsin{x} &= \int \frac{1}{\sqrt{1 - x^{2}}} \mathrm{d} x \\ \frac{\mathrm{d}}{\mathrm{d} x} \arccos{x} &= - \frac{1}{\sqrt{1 - x^{2}}} \\ \arccos{x} &= - \int \frac{1}{\sqrt{1 - x^{2}}} \mathrm{d} x \\ \frac{\mathrm{d}}{\mathrm{d} x} \arctan{x} &= \frac{1}{1 + x^{2}} \\ \arctan{x} &= \int \frac{1}{1 + x^{2}} \mathrm{d} x \end{aligned} \]

那么代入 \(f\left(x\right)\) 就有:

\[ \begin{aligned} \frac{\mathrm{d}}{\mathrm{d} x} \arcsin{f\left(x\right)} &= \frac{f'\left(x\right)}{\sqrt{1 - f^{2}\left(x\right)}} \\ \arcsin{f\left(x\right)} &= \int \frac{f'\left(x\right)}{\sqrt{1 - f^{2}\left(x\right)}} \mathrm{d} x \\ \frac{\mathrm{d}}{\mathrm{d} x} \arccos{f\left(x\right)} &= - \frac{f'\left(x\right)}{\sqrt{1 - f^{2}\left(x\right)}} \\ \arccos{f\left(x\right)} &= - \int \frac{f'\left(x\right)}{\sqrt{1 - f^{2}\left(x\right)}} \mathrm{d} x \\ \frac{\mathrm{d}}{\mathrm{d} x} \arctan{f\left(x\right)} &= \frac{f'\left(x\right)}{1 + f^{2}\left(x\right)} \\ \arctan{f\left(x\right)} &= \int \frac{f'\left(x\right)}{1 + f^{2}\left(x\right)} \mathrm{d} x \end{aligned} \]

直接按式子求就可以了。(我也不清楚Luogu凭什么给这个玩意评黑)

实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
void Asin(int n, ll *F, ll *G)
{
    static ll A[MAXN], B[MAXN];
    copy(F, F + n, A);
    int lim;
    for (lim = 1; lim < (n << 1); lim <<= 1);
    for (int i = 0; i < lim; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (lim >> 1));
    NTT(A, lim, 1);
    for (int i = 0; i < lim; i++) A[i] = A[i] * A[i] % p;
    NTT(A, lim, -1), fill(A + n, A + lim, 0);
    for (int i = 0; i < n; i++) A[i] = (p - A[i]) % p;
    A[0] = (A[0] + 1) % p, Sqrt(n, A, B), fill(A, A + n, 0), Inv(n, B, A),
    fill(B, B + n, 0), Dert(n, F, B);
    for (int i = 0; i < lim; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (lim >> 1));
    NTT(A, lim, 1), NTT(B, lim, 1);
    for (int i = 0; i < lim; i++) A[i] = A[i] * B[i] % p;
    NTT(A, lim, -1), fill(A + n, A + lim, 0), fill(B, B + lim, 0), Int(n, A, G),
        fill(A, A + n, 0);
}
void Atan(int n, ll *F, ll *G)
{
    static ll A[MAXN], B[MAXN];
    copy(F, F + n, A);
    int lim;
    for (lim = 1; lim < (n << 1); lim <<= 1);
    for (int i = 0; i < lim; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (lim >> 1));
    NTT(A, lim, 1);
    for (int i = 0; i < lim; i++) A[i] = A[i] * A[i] % p;
    NTT(A, lim, -1), fill(A + n, A + lim, 0),
        A[0] = (A[0] + 1) % p, Inv(n, A, B), fill(A, A + n, 0), Dert(n, F, A);
    for (int i = 0; i < lim; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (lim >> 1));
    NTT(A, lim, 1), NTT(B, lim, 1);
    for (int i = 0; i < lim; i++) A[i] = A[i] * B[i] % p;
    NTT(A, lim, -1), fill(A + n, A + lim, 0), fill(B, B + lim, 0), Int(n, A, G),
        fill(A, A + n, 0);
}

总结

下面的代码实现了 Poly 类,实现了本页面的所有操作,但常数很大。

实现
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
#include <bits/stdc++.h>
using namespace std;
#define endl '\n'
typedef long long ll;
const int MAXN = 6e5 + 10;
const ll p = 998244353, g = 3;
inline ll qpow(ll a, ll b)
{
    static ll res;
    for (res = 1; b; b >>= 1, a = a * a % p)
        if (b & 1) res = res * a % p;
    return res;
}
inline ll getinv(const ll x) { return qpow(x, p - 2); }
const ll gi = getinv(g), img = qpow(g, (p - 1) / 4), inv2 = getinv(2), inv2img = getinv(img * 2);
int rev[MAXN];
ll inv[MAXN];
inline void initinv(const int n)
{
    inv[1] = 1;
    for (int i = 2; i <= n; i++) inv[i] = (p - p / i) * inv[p % i] % p;
}
inline int getrev(const int n)
{
    int lim = 1;
    while (lim <= n) lim <<= 1;
    for (int i = 0; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (lim >> 1));
    return lim;
}
struct Poly : public vector<ll>
{
    using vector<ll>::vector; // 继承 vector<ll> 的所有构造函数
    static inline void read(const int n, Poly &F)
    {
        F.resize(n);
        for (ll &x : F) cin >> x, x %= p;
    }
    static inline void write(const Poly F)
    {
        for (ll x : F) cout << x << ' ';
        cout << endl;
    }
    friend inline Poly operator+(Poly F, Poly G)
    {
        if (F.size() > G.size()) std::swap(F, G);
        for (int i = 0; i < F.size(); i++) G[i] = (G[i] + F[i]) % p;
        return G;
    }
    friend inline Poly operator-(Poly F, const Poly G)
    {
        if (F.size() < G.size()) F.resize(G.size());
        for (int i = 0; i < G.size(); i++) F[i] = (F[i] - G[i] + p) % p;
        return F;
    }
    friend inline Poly operator*(Poly F, ll x)
    {
        x = (x + p) % p;
        for (ll &y : F) y = x * y % p;
        return F;
    }
    friend inline Poly operator*(const ll x, const Poly F) { return F * x; }
    static inline Poly NTT(const Poly F, const int n, const int op)
    {
        Poly A(F);
        A.resize(n);
        for (int i = 0; i < n; i++)
            if (i < rev[i]) std::swap(A[i], A[rev[i]]);
        for (int i = 2; i <= n; i <<= 1)
        {
            ll g1 = qpow(op == 1 ? g : gi, p / i);
            for (int j = 0; j < n; j += i)
            {
                ll gk = 1;
                for (int k = j; k < j + i / 2; k++)
                {
                    ll x = A[k], y = A[k + i / 2] * gk % p;
                    A[k] = (x + y) % p;
                    A[k + i / 2] = (x - y + p) % p;
                    gk = gk * g1 % p;
                }
            }
        }
        if (op == 1) return A;
        const ll ni = getinv(n);
        for (int i = 0; i < n; i++) A[i] = A[i] * ni % p;
        return A;
    }
    friend inline Poly operator*(Poly F, Poly G)
    {
        const int n = F.size(), m = G.size(), len = n + m - 1, lim = getrev(len);
        F = NTT(F, lim, 1), G = NTT(G, lim, 1);
        for (int i = 0; i < lim; i++) F[i] = F[i] * G[i] % p;
        return F = NTT(F, lim, -1), F.resize(len), F;
    }
    static inline Poly Dert(const Poly F)
    {
        const int n = F.size();
        Poly G(n - 1);
        for (int i = 1; i < n; i++) G[i - 1] = F[i] * i % p;
        return G;
    }
    static inline Poly Int(const Poly F)
    {
        const int n = F.size();
        Poly G(n + 1);
        for (int i = 0; i < n; i++) G[i + 1] = F[i] * inv[i + 1] % p;
        return G;
    }
    static inline Poly Inv(Poly F)
    {
        int n = F.size(), len, lim = getrev(n);
        Poly G{getinv(F[0])};
        F.resize(lim);
        for (len = 1; len < (n << 1); len <<= 1)
        {
            G.resize(lim = getrev(len << 1));
            Poly A(F.begin(), F.begin() + len), B(G.begin(), G.begin() + len);
            A = NTT(A, lim, 1), B = NTT(B, lim, 1);
            for (int i = 0; i < lim; i++) G[i] = ((2 - B[i] * A[i] % p) % p * B[i] % p + p) % p;
            G = NTT(G, lim, -1), G.resize(len);
        }
        return G.resize(n), G;
    }
    static inline void Div(Poly F, Poly G, Poly &Q, Poly &R)
    {
        const int n = F.size(), m = G.size(), len = n - m + 1;
        Poly A(F), B(G);
        reverse(A.begin(), A.end()), reverse(B.begin(), B.end()), A.resize(len), B.resize(len), Q = A * Inv(B), Q.resize(len), reverse(Q.begin(), Q.end()), R = F - G * Q, R.resize(m - 1);
    }
    static inline Poly Ln(const Poly F)
    {
        Poly G = Dert(F) * Inv(F);
        return G.resize(F.size() - 1), Int(G);
    }
    static inline Poly Exp(Poly F)
    {
        int n = F.size(), len, lim = getrev(n);
        Poly G{1};
        F.resize(lim);
        for (len = 1; len < (n << 1); len <<= 1)
        {
            G.resize(lim = getrev(len << 1));
            Poly A(F.begin(), F.begin() + len), B(G.begin(), G.begin() + len);
            G = B * (Poly({1}) - Ln(B) + A), G.resize(len);
        }
        return G.resize(n), G;
    }
    static inline Poly Sqrt(Poly F)
    {
        int n = F.size(), len, lim = getrev(n);
        Poly G{1};
        F.resize(lim);
        for (len = 1; len < (n << 1); len <<= 1)
        {
            G.resize(lim = getrev(len << 1));
            Poly A(F.begin(), F.begin() + len), B(G.begin(), G.begin() + len), I = Inv(B);
            A = NTT(A, lim, 1), B = NTT(B, lim, 1), I = NTT(I, lim, 1);
            for (int i = 0; i < lim; i++) G[i] = (B[i] * B[i] % p + A[i]) % p * inv2 % p * I[i] % p;
            G = NTT(G, lim, -1), G.resize(len);
        }
        return G.resize(n), G;
    }
    static inline Poly Pow(const Poly F, const ll k) { return Exp(Ln(F) * k); }
    static inline Poly Sin(const Poly F) { return (Exp(img * F) - Exp((p - img) * F)) * inv2img; }
    static inline Poly Cos(const Poly F) { return (Exp(img * F) + Exp((p - img) * F)) * inv2; }
    static inline Poly Asin(const Poly F)
    {
        Poly G = Poly({1}) - F * F;
        return G.resize(F.size()), G = (Int(Dert(F) * Inv(Sqrt(G)))), G.resize(F.size()), G;
    }
    static inline Poly Atan(const Poly F)
    {
        Poly G = Poly({1}) + F * F;
        return G.resize(F.size()), G = (Int(Dert(F) * Inv(G))), G.resize(F.size()), G;
    }
};
int n, m;
int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    // Usage:
    initinv(6e5);
    cin >> n >> m;
    Poly F, G, Q, R;
    Poly::read(n, F), Poly::read(m, G);
    Poly::write(F * G);
    Poly::write(Poly::Inv(F));
    Poly::write(Poly::Sqrt(F));
    Poly::Div(F, G, Q, R), Poly::write(Q), Poly::write(R);
    Poly::write(Poly::Ln(F));
    Poly::write(Poly::Exp(F));
    Poly::write(Poly::Sin(F));
    Poly::write(Poly::Cos(F));
    Poly::write(Poly::Asin(F));
    Poly::write(Poly::Atan(F));
    return 0;
}