Elementary Polynomial Functions

This page covers common elementary operations on polynomials. Specifically, it includes the following:

  1. Polynomial inversion
  2. Polynomial square roots
  3. Polynomial division
  4. Polynomial remainder
  5. Polynomial exponentials
  6. Polynomial logarithms
  7. Polynomial trigonometric functions
  8. Polynomial inverse trigonometric functions
Elementary and non-elementary functions

The definition of an elementary function is as follows1:

If there exists a mapping u→∂uu\to \partial u in a field FF satisfying:

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

then this field is called a differential field.

If a function uu on the differential field FF satisfies any one of the following conditions, then uu is called an elementary function:

  1. uu is an algebraic function over FF.
  2. uu is an exponential function over FF, meaning that there exists a∈Fa\in F such that ∂u=u∂a\partial u=u\partial a.
  3. uu is a logarithmic function over FF, meaning that there exists a∈Fa\in F such that ∂u=∂aa\partial u=\frac{\partial a}{a}.

The following are common elementary functions:

  1. Algebraic functions: functions f(x)f(x) for which there exists a finite-degree polynomial PP such that P(f(x))=0P(f(x))=0, such as 2x+12x+1,x\sqrt{x},(1+x2)−1(1+x^2)^{-1},∣x∣|x|.

  2. Exponential functions

  3. Logarithmic functions

  4. Trigonometric functions

  5. Inverse trigonometric functions

  6. Hyperbolic functions

  7. Inverse hyperbolic functions

  8. Compositions of the above functions, such as:

    etan⁡x1+x2sin⁡(1+ln⁡2x)\frac{\mathrm{e}^{\tan x}}{1+x^2}\sin\left(\sqrt{1+\ln^2 x}\right) −iln⁡(x+i1−x2)-\mathrm{i} \ln\left(x+\mathrm{i}\sqrt{1-x^2}\right)

The following is a common non-elementary function:

  1. The error function:

    erf⁡(x):=2π∫0xexp⁡(−t2)dt\operatorname{erf}(x):=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp\left(-t^2\right)\mathrm{d}t

Polynomial Inversion

Given a polynomial f(x)f\left(x\right), find f−1(x)f^{-1}\left(x\right).

Solution

Doubling Method

First, it is easy to see that

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

Assume that the inverse f0−1(x)f^{-1}_{0}\left(x\right) of f(x)f\left(x\right) modulo x⌈n2⌉x^{\left\lceil\frac{n}{2}\right\rceil} has already been found. Then:

f(x)f0−1(x)≡1(modx⌈n2⌉)f(x)f−1(x)≡1(modx⌈n2⌉)f−1(x)−f0−1(x)≡0(modx⌈n2⌉)\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}

Squaring both sides gives:

f−2(x)−2f−1(x)f0−1(x)+f0−2(x)≡0(modxn)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}}

Multiplying both sides by f(x)f\left(x\right) and rearranging gives:

f−1(x)≡f0−1(x)(2−f(x)f0−1(x))(modxn)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}}

The result follows by recursion.

Time complexity

T(n)=T(n2)+O(nlog⁡n)=O(nlog⁡n)T\left(n\right)=T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log{n}\right)

Newton's Method

See Newton's Method.

Graeffe's Method

To find f−1(x) mod x2nf^{-1}(x)\bmod x^{2n}, consider

f−1(x) mod x2n=f(−x)(f(x)f(−x))−1 mod x2n=f(−x)g−1(x2) mod x2n\begin{aligned} f^{-1}(x)\bmod x^{2n}&= f(-x)(f(x)f(-x))^{-1}\bmod x^{2n}\\ &=f(-x)g^{-1}(x^2)\bmod x^{2n} \end{aligned}

It is enough to find g−1(x) mod xng^{-1}(x)\bmod x^n to recover g−1(x2) mod x2ng^{-1}(x^2)\bmod x^{2n}, because f(x)f(−x)f(x)f(-x) is an even function. The time complexity is the same as above.

Code

Polynomial Inversion
constexpr int MAXN = 262144;
constexpr int mod = 998244353;

using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;

void polyinv(const poly &h, const int n, poly &f) {
  /* f = 1 / h = f_0 (2 - f_0 h) */
  static poly_t inv_t;
  std::fill(f, f + n + n, 0);
  f[0] = fpow(h[0], mod - 2);
  for (int t = 2; t <= n; t <<= 1) {
    const int t2 = t << 1;
    std::copy(h, h + t, inv_t);
    std::fill(inv_t + t, inv_t + t2, 0);

    DFT(f, t2);
    DFT(inv_t, t2);
    for (int i = 0; i != t2; ++i)
      f[i] = (i64)f[i] * sub(2, (i64)f[i] * inv_t[i] % mod) % mod;
    IDFT(f, t2);

    std::fill(f + t, f + t2, 0);
  }
}

Example

  1. Counting labeled simple undirected connected graphs: 「POJ 1737」Connected Graph

Polynomial Square Roots

Given a polynomial g(x)g\left(x\right), find f(x)f\left(x\right) satisfying:

f2(x)≡g(x)(modxn)f^{2}\left(x\right)\equiv g\left(x\right) \pmod{x^{n}}

Solution

Doubling Method

First consider the case where [x0]g(x)\left[x^0\right]g(x) is not 00.

It is easy to see that:

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

If [x0]g(x)\left[x^0\right]g(x) has no square root, then the polynomial g(x)g(x) has no square root.

[x0]g(x)\left[x^0\right]g(x) may have multiple square roots; choosing different roots yields different f(x)f(x).

Assume that the square root f0(x)f_{0}\left(x\right) of g(x)g\left(x\right) modulo x⌈n2⌉x^{\left\lceil\frac{n}{2}\right\rceil} has already been found. Then:

f02(x)≡g(x)(modx⌈n2⌉)f02(x)−g(x)≡0(modx⌈n2⌉)(f02(x)−g(x))2≡0(modxn)(f02(x)+g(x))2≡4f02(x)g(x)(modxn)(f02(x)+g(x)2f0(x))2≡g(x)(modxn)f02(x)+g(x)2f0(x)≡f(x)(modxn)2−1f0(x)+2−1f0−1(x)g(x)≡f(x)(modxn)\begin{aligned} f_{0}^{2}\left(x\right)&\equiv g\left(x\right) &\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\\ f_{0}^{2}\left(x\right)-g\left(x\right)&\equiv 0 &\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\\ \left(f_{0}^{2}\left(x\right)-g\left(x\right)\right)^{2}&\equiv 0 &\pmod{x^{n}}\\ \left(f_{0}^{2}\left(x\right)+g\left(x\right)\right)^{2}&\equiv 4f_{0}^{2}\left(x\right)g\left(x\right) &\pmod{x^{n}}\\ \left(\frac{f_{0}^{2}\left(x\right)+g\left(x\right)}{2f_{0}\left(x\right)}\right)^{2}&\equiv g\left(x\right) &\pmod{x^{n}}\\ \frac{f_{0}^{2}\left(x\right)+g\left(x\right)}{2f_{0}\left(x\right)}&\equiv f\left(x\right) &\pmod{x^{n}}\\ 2^{-1}f_{0}\left(x\right)+2^{-1}f_{0}^{-1}\left(x\right)g\left(x\right)&\equiv f\left(x\right) &\pmod{x^{n}} \end{aligned}

Compute by doubling.

Time complexity

T(n)=T(n2)+O(nlog⁡n)=O(nlog⁡n)T\left(n\right)=T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log{n}\right)

Another version with a smaller constant is to maintain f−1(x)f^{-1}\left(x\right) alongside f(x)f\left(x\right) during doubling, instead of computing the inverse each time.

When [x0]g(x)≠1\left[x^{0}\right]g\left(x\right)\neq 1, a quadratic-residue algorithm may be needed to compute [x0]f(x)\left[x^{0}\right]f\left(x\right).

The method above requires knowing the inverse of f0(x)f_{0}(x), so its constant term cannot be 00.

If [x0]g(x)=0\left[x^0\right]g(x) = 0, factor g(x)g(x) as xkh(x)x^{k}h(x), where [x0]h(x)≠0\left[x^0\right]h(x) \not = 0.

  • If kk is odd, then g(x)g(x) has no square root.

  • If kk is even, find the square root h(x)\sqrt{h(x)} of h(x)h(x), and then obtain f(x)≡xk/2h(x)(modxn)f(x) \equiv x^{k/2} \sqrt{h(x)} \pmod{x^{n}}.

Luogu Template Problem P5205 【Template】Polynomial Square Root Reference Code
#include <algorithm>
#include <iostream>
using namespace std;

constexpr int MAXN = 1 << 20, mod = 998244353;

int a[MAXN], b[MAXN], g[MAXN], gg[MAXN];

int qpow(int x, int y) {  // fast exponentiation
  int ans = 1;

  while (y) {
    if (y & 1) {
      ans = (long long)1 * ans * x % mod;
    }
    x = (long long)1 * x * x % mod;
    y >>= 1;
  }
  return ans;
}

int inv2 = qpow(2, mod - 2);  // modular inverse

void change(int *f, int len) {
  for (int i = 1, j = len / 2; i < len - 1; i++) {
    if (i < j) {
      swap(f[i], f[j]);
    }

    int k = len / 2;
    while (j >= k) {
      j -= k;
      k /= 2;
    }
    if (j < k) {
      j += k;
    }
  }
}

void NTT(int *f, int len, int type) {  // NTT
  change(f, len);

  for (int q = 2; q <= len; q <<= 1) {
    int nxt = qpow(3, (mod - 1) / q);
    for (int i = 0; i < len; i += q) {
      int w = 1;

      for (int k = i; k < i + (q >> 1); k++) {
        int x = f[k];
        int y = (long long)1 * w * f[k + (q / 2)] % mod;

        f[k] = (x + y) % mod;
        f[k + (q / 2)] = (x - y + mod) % mod;
        w = (long long)1 * w * nxt % mod;
      }
    }
  }

  if (type == -1) {
    reverse(f + 1, f + len);
    int iv = qpow(len, mod - 2);

    for (int i = 0; i < len; i++) {
      f[i] = (long long)1 * f[i] * iv % mod;
    }
  }
}

void inv(int deg, int *f, int *h) {  // polynomial inversion
  if (deg == 1) {
    h[0] = qpow(f[0], mod - 2);
    return;
  }

  inv((deg + 1) >> 1, f, h);

  int len = 1;
  while (len < deg * 2) {  // doubling
    len *= 2;
  }

  copy(f, f + deg, gg);
  fill(gg + deg, gg + len, 0);

  NTT(gg, len, 1);
  NTT(h, len, 1);
  for (int i = 0; i < len; i++) {
    h[i] = (long long)1 * (2 - (long long)1 * gg[i] * h[i] % mod + mod) % mod *
           h[i] % mod;
  }

  NTT(h, len, -1);
  fill(h + deg, h + len, 0);
}

int n, t[MAXN];

// deg: degree
// f: input array whose square root is taken
// h: output array
void sqrt(int deg, int *f, int *h) {
  if (deg == 1) {
    h[0] = 1;
    return;
  }

  sqrt((deg + 1) >> 1, f, h);

  int len = 1;
  while (len < deg * 2) {  // doubling
    len *= 2;
  }
  fill(g, g + len, 0);
  inv(deg, h, g);
  copy(f, f + deg, t);
  fill(t + deg, t + len, 0);
  NTT(t, len, 1);
  NTT(g, len, 1);
  NTT(h, len, 1);

  for (int i = 0; i < len; i++) {
    h[i] = (long long)1 * inv2 *
           ((long long)1 * h[i] % mod + (long long)1 * g[i] * t[i] % mod) % mod;
  }
  NTT(h, len, -1);
  fill(h + deg, h + len, 0);
}

int main() {
  cin >> n;

  for (int i = 0; i < n; i++) {
    cin >> a[i];
  }
  sqrt(n, a, b);

  for (int i = 0; i < n; i++) {
    cout << b[i] << ' ';
  }

  return 0;
}

Newton's Method

See Newton's Method.

Example

  1. 「Codeforces Round #250」E. The Child and Binary Tree

Polynomial Division & Modulo

Given polynomials f(x),g(x)f\left(x\right),g\left(x\right), find the quotient Q(x)Q\left(x\right) and remainder R(x)R\left(x\right) when f(x)f\left(x\right) is divided by g(x)g\left(x\right).

Solution

If the effect of R(x)R\left(x\right) can be eliminated, the problem can be solved directly using polynomial inversion.

Consider constructing the transformation

fR(x)=xdeg⁡ff(1x)f^{R}\left(x\right)=x^{\operatorname{deg}{f}}f\left(\frac{1}{x}\right)

It can be seen that this simply reverses the coefficients of f(x)f\left(x\right).

Let n=deg⁡f,m=deg⁡gn=\operatorname{deg}{f},m=\operatorname{deg}{g}.

Substitute 1x\frac{1}{x} for xx in f(x)=Q(x)g(x)+R(x)f\left(x\right)=Q\left(x\right)g\left(x\right)+R\left(x\right) and multiply both sides by xnx^{n} to obtain:

xnf(1x)=xn−mQ(1x)xmg(1x)+xn−m+1xm−1R(1x)fR(x)=QR(x)gR(x)+xn−m+1RR(x)\begin{aligned} 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}\left(x\right)&=Q^{R}\left(x\right)g^{R}\left(x\right)+x^{n-m+1}R^{R}\left(x\right) \end{aligned}

Notice that the coefficients of RR(x)R^{R}\left(x\right) in the equation above are multiplied by xn−m+1x^{n-m+1}, so taking the equation modulo xn−m+1x^{n-m+1} eliminates the effect of RR(x)R^{R}\left(x\right).

Also, since the degree of QR(x)Q^{R}\left(x\right) is (n−m)<(n−m+1)\left(n-m\right)<\left(n-m+1\right), QR(x)Q^{R}\left(x\right) is unaffected.

Therefore:

fR(x)≡QR(x)gR(x)(modxn−m+1)f^{R}\left(x\right)\equiv Q^{R}\left(x\right)g^{R}\left(x\right)\pmod{x^{n-m+1}}

The quotient Q(x)Q\left(x\right) can be obtained using polynomial inversion, and substituting it back gives R(x)R\left(x\right).

Time complexity O(nlog⁡n)O\left(n\log{n}\right).

Polynomial Logarithms & Exponentials

Given a polynomial f(x)f(x), find ln⁡f(x)\ln{f(x)} and exp⁡f(x)\exp{f(x)} modulo xnx^{n}.

Solution

Ordinary Method

Polynomial logarithm

First, for a polynomial f(x)f(x), if ln⁡f(x)\ln{f(x)} exists, then by its definition, it must satisfy:

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

Differentiate ln⁡f(x)\ln{f(x)} and then integrate to obtain:

dln⁡f(x)dx≡f′(x)f(x)(modxn)ln⁡f(x)≡∫dln⁡f(x)≡∫f′(x)f(x)dx(modxn)\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 \mathrm{d} \ln{f(x)} \equiv \int\frac{f'(x)}{f(x)} \mathrm{d} x & \pmod{x^{n}} \end{aligned}

Polynomial differentiation and integration each take O(n)O(n) time, while inversion takes O(nlog⁡n)O(n\log{n}) time. Thus, computing the polynomial ln⁡\ln takes O(nlog⁡n)O(n\log{n}) time.

Polynomial exponential

First, for a polynomial f(x)f(x), if exp⁡f(x)\exp{f(x)} exists, then it must satisfy:

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

Otherwise, the constant term of exp⁡f(x)\exp{f(x)} does not converge.

Differentiate exp⁡f(x)\exp{f(x)} to obtain:

dexp⁡f(x)dx≡exp⁡f(x)f′(x)(modxn)\frac{\mathrm{d} \exp{f(x)}}{\mathrm{d} x} \equiv \exp{f(x)}f'(x)\pmod{x^{n}}

Comparing coefficients on both sides gives:

[xn−1]dexp⁡f(x)dx=∑i=0n−1([xi]exp⁡f(x))([xn−i−1]f′(x))[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[xn]exp⁡f(x)=∑i=0n−1([xi]exp⁡f(x))((n−i)[xn−i]f(x))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)

This can be solved using divide-and-conquer FFT.

Time complexity O(nlog⁡2n)O(n\log^{2}{n}).

Newton's Method

Using Newton's Method, the polynomial exp⁡\exp can be computed in O(nlog⁡n)O(n\log{n}) time.

Code

Polynomial ln/exp
constexpr int MAXN = 262144;
constexpr int mod = 998244353;

using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;

void derivative(const poly &h, const int n, poly &f) {
  for (int i = 1; i != n; ++i) f[i - 1] = (i64)h[i] * i % mod;
  f[n - 1] = 0;
}

void integrate(const poly &h, const int n, poly &f) {
  for (int i = n - 1; i; --i) f[i] = (i64)h[i - 1] * inv[i] % mod;
  f[0] = 0; /* C */
}

void polyln(const poly &h, const int n, poly &f) {
  /* f = ln h = ∫ h' / h dx */
  assert(h[0] == 1);
  static poly_t ln_t;
  const int t = n << 1;

  derivative(h, n, ln_t);
  std::fill(ln_t + n, ln_t + t, 0);
  polyinv(h, n, f);

  DFT(ln_t, t);
  DFT(f, t);
  for (int i = 0; i != t; ++i) ln_t[i] = (i64)ln_t[i] * f[i] % mod;
  IDFT(ln_t, t);

  integrate(ln_t, n, f);
}

void polyexp(const poly &h, const int n, poly &f) {
  /* f = exp(h) = f_0 (1 - ln f_0 + h) */
  assert(h[0] == 0);
  static poly_t exp_t;
  std::fill(f, f + n + n, 0);
  f[0] = 1;
  for (int t = 2; t <= n; t <<= 1) {
    const int t2 = t << 1;

    polyln(f, t, exp_t);
    exp_t[0] = sub(pls(h[0], 1), exp_t[0]);
    for (int i = 1; i != t; ++i) exp_t[i] = sub(h[i], exp_t[i]);
    std::fill(exp_t + t, exp_t + t2, 0);

    DFT(f, t2);
    DFT(exp_t, t2);
    for (int i = 0; i != t2; ++i) f[i] = (i64)f[i] * exp_t[i] % mod;
    IDFT(f, t2);

    std::fill(f + t, f + t2, 0);
  }
}

Example

  1. Computing fk(x)f^{k}(x)

    The ordinary approach is polynomial fast exponentiation, with time complexity O(nlog⁡nlog⁡k)O(n\log{n}\log{k}).

    When [x0]f(x)=1[x^{0}]f(x)=1, we have:

    fk(x)=exp⁡(kln⁡f(x))f^{k}(x)=\exp{\left(k\ln{f(x)}\right)}

    When [x0]f(x)≠1[x^{0}]f(x)\neq 1, let the lowest-degree term of f(x)f(x) be fixif_{i}x^{i}. Then:

    fk(x)=fikxikexp⁡(kln⁡f(x)fixi)f^{k}(x)=f_{i}^{k}x^{ik}\exp{\left(k\ln{\frac{f(x)}{f_{i}x^{i}}}\right)}

    Time complexity O(nlog⁡n)O(n\log{n}).

Polynomial Trigonometric Functions

Given a polynomial f(x)f\left(x\right), find sin⁡f(x),cos⁡f(x)\sin{f\left(x\right)}, \cos{f\left(x\right)}, and tan⁡f(x)\tan{f\left(x\right)} modulo xnx^{n}.

Solution

By Euler's formula, (eix=cos⁡x+isin⁡x)\left(\mathrm{e}^{\mathrm{i}x} = \cos{x} + \mathrm{i}\sin{x}\right), we obtain another expression for trigonometric functions:

sin⁡x=eix−e−ix2icos⁡x=eix+e−ix2\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}

Substituting f(x)f\left(x\right) gives:

sin⁡f(x)=exp⁡(if(x))−exp⁡(−if(x))2icos⁡f(x)=exp⁡(if(x))+exp⁡(−if(x))2\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}

Writing a program directly from the expressions above gives sin⁡f(x)\sin{f\left(x\right)} and cos⁡f(x)\cos{f\left(x\right)} modulo xnx^{n}. Then tan⁡f(x)\tan{f\left(x\right)} can be obtained from tan⁡f(x)=sin⁡f(x)cos⁡f(x)\tan{f\left(x\right)} = \frac{\sin{f\left(x\right)}}{\cos{f\left(x\right)}}.

Code

Polynomial Trigonometric Functions

Notice that we are performing NTT over Z998244353\mathbb{Z}_{998244353}. Accordingly, the imaginary unit i\mathrm{i} should be replaced by 8658371886583718 or 911660635911660635:

i=−1≡998244352(mod998244353)  ⟹  ori≡86583718(mod998244353)ori≡911660635(mod998244353)\begin{aligned} & \mathrm{i} = \sqrt{-1} \equiv \sqrt{998244352} \pmod{998244353} \\ \implies & \phantom{\text{or}} \quad \mathrm{i} \equiv 86583718 \pmod{998244353} \\ & \text{or} \quad \mathrm{i} \equiv 911660635 \pmod{998244353} \end{aligned}
constexpr int MAXN = 262144;
constexpr int mod = 998244353;
constexpr int imgunit = 86583718; /* sqrt(-1) = sqrt(998233452) */

using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;

void polytri(const poly &h, const int n, poly &sin_t, poly &cos_t) {
  /* sin(f) = (exp(i * f) - exp(- i * f)) / 2i */
  /* cos(f) = (exp(i * f) + exp(- i * f)) / 2 */
  /* tan(f) = sin(f) / cos(f) */
  assert(h[0] == 0);
  static poly_t tri1_t, tri2_t;

  for (int i = 0; i != n; ++i) tri2_t[i] = (i64)h[i] * imgunit % mod;
  polyexp(tri2_t, n, tri1_t);
  polyinv(tri1_t, n, tri2_t);

  if (sin_t != nullptr) {
    const int invi = fpow(pls(imgunit, imgunit), mod - 2);
    for (int i = 0; i != n; ++i)
      sin_t[i] = (i64)(tri1_t[i] - tri2_t[i] + mod) * invi % mod;
  }
  if (cos_t != nullptr) {
    for (int i = 0; i != n; ++i) cos_t[i] = div2(pls(tri1_t[i], tri2_t[i]));
  }
}

Polynomial Inverse Trigonometric Functions

Given a polynomial f(x)f\left(x\right), find arcsin⁡f(x),arccos⁡f(x)\arcsin{f\left(x\right)}, \arccos{f\left(x\right)}, and arctan⁡f(x)\arctan{f\left(x\right)} modulo xnx^{n}.

Solution

Following the method for computing the polynomial ln⁡\ln, differentiate and then integrate the inverse trigonometric functions to obtain:

ddxarcsin⁡x=11−x2arcsin⁡x=∫11−x2dxddxarccos⁡x=−11−x2arccos⁡x=−∫11−x2dxddxarctan⁡x=11+x2arctan⁡x=∫11+x2dx\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}

Substituting f(x)f\left(x\right) gives:

ddxarcsin⁡f(x)=f′(x)1−f2(x)arcsin⁡f(x)=∫f′(x)1−f2(x)dxddxarccos⁡f(x)=−f′(x)1−f2(x)arccos⁡f(x)=−∫f′(x)1−f2(x)dxddxarctan⁡f(x)=f′(x)1+f2(x)arctan⁡f(x)=∫f′(x)1+f2(x)dx\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}

The result follows by evaluating these expressions directly.

Code

Polynomial Inverse Trigonometric Functions
constexpr int MAXN = 262144;
constexpr int mod = 998244353;

using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;

void derivative(const poly &h, const int n, poly &f) {
  for (int i = 1; i != n; ++i) f[i - 1] = (i64)h[i] * i % mod;
  f[n - 1] = 0;
}

void integrate(const poly &h, const int n, poly &f) {
  for (int i = n - 1; i; --i) f[i] = (i64)h[i - 1] * inv[i] % mod;
  f[0] = 0; /* C */
}

void polyarcsin(const poly &h, const int n, poly &f) {
  /* arcsin(f) = ∫ f' / sqrt(1 - f^2) dx  */
  static poly_t arcsin_t;
  const int t = n << 1;
  std::copy(h, h + n, arcsin_t);
  std::fill(arcsin_t + n, arcsin_t + t, 0);

  DFT(arcsin_t, t);
  for (int i = 0; i != t; ++i) arcsin_t[i] = sqr(arcsin_t[i]);
  IDFT(arcsin_t, t);

  arcsin_t[0] = sub(1, arcsin_t[0]);
  for (int i = 1; i != n; ++i)
    arcsin_t[i] = arcsin_t[i] ? mod - arcsin_t[i] : 0;

  polysqrt(arcsin_t, n, f);
  polyinv(f, n, arcsin_t);
  derivative(h, n, f);

  DFT(f, t);
  DFT(arcsin_t, t);
  for (int i = 0; i != t; ++i) arcsin_t[i] = (i64)f[i] * arcsin_t[i] % mod;
  IDFT(arcsin_t, t);

  integrate(arcsin_t, n, f);
}

void polyarccos(const poly &h, const int n, poly &f) {
  /* arccos(f) = - ∫ f' / sqrt(1 - f^2) dx  */
  polyarcsin(h, n, f);
  for (int i = 0; i != n; ++i) f[i] = f[i] ? mod - f[i] : 0;
}

void polyarctan(const poly &h, const int n, poly &f) {
  /* arctan(f) = ∫ f' / (1 + f^2) dx  */
  static poly_t arctan_t;
  const int t = n << 1;
  std::copy(h, h + n, arctan_t);
  std::fill(arctan_t + n, arctan_t + t, 0);

  DFT(arctan_t, t);
  for (int i = 0; i != t; ++i) arctan_t[i] = sqr(arctan_t[i]);
  IDFT(arctan_t, t);

  inc(arctan_t[0], 1);
  std::fill(arctan_t + n, arctan_t + t, 0);

  polyinv(arctan_t, n, f);
  derivative(h, n, arctan_t);

  DFT(f, t);
  DFT(arctan_t, t);
  for (int i = 0; i != t; ++i) arctan_t[i] = (i64)f[i] * arctan_t[i] % mod;
  IDFT(arctan_t, t);

  integrate(arctan_t, n, f);
}

Footnotes

  1. Elementary function——Wikipedia ↩

This content comes from OI Wiki: Elementary Polynomial Functions. See the original for authors and revision history. AIOJ has adapted the formatting, in-site links, and resource paths.

Content is shared under CC BY-SA 4.0 and SATA license; additional terms are governed by the original notice. Thanks to the OI Wiki authors and contributors; per the SATA terms, please visit the OI Wiki repository and star it to show your support.