More on D-Finite Functions I: Theories

#	User	Rating
1	jiangly	3977
2	tourist	3815
3	jqdai0815	3682
4	ksun48	3614
5	orzdevinwang	3526
6	ecnerwala	3514
7	Benq	3483
8	hos.lyric	3381
9	gamegame	3374
10	heuristica	3358

#	User	Contrib.
1	cry	170
2	-is-this-fft-	162
2	Um_nik	162
4	atcoder_official	160
5	djm03178	157
5	Dominater069	157
7	adamant	154
8	luogu_official	152
8	awoo	152
10	TheScrasse	147

Hello everyone.

Several years ago jqdai0815 released a blog on the topic of D-Finite functions. In this blog, I will introduce some more theories, and show the relation between some CP problems. Also, I found out that what I show is like opening Pandora's box: it provides some new ways to solve problems, but they also somehow prevent us from really understanding them. Therefore, I'll also talk about some of my personal critiques on this topic.

Since this blog is quite long, I'll divide it into several parts. This is the first part.

Appetizer

We start with this toy problem.

$$$n$$$-King Problem. Given $$$n$$$, count the number of permutations $$$\sigma$$$ of $$$[n]$$$ such that for any $$$i$$$, $$$|\sigma(i) - \sigma(i+1)| \neq 1$$$.

Of course, this sequence has already been studied before, it is the A002464 in OEIS. You can find its recurrence relation

$$$ a_n = (n+1)a_{n-1} - (n-2)a_{n-2} - (n-5)a_{n-3}+ (n-3)a_{n-4}. $$$

But how to prove this? A combinatorial proof is possible, but not that easy. Here I will show a different way to prove this.

We first consider the inclusion-exclusion principle. Suppose we force a subset $$$S$$$ of $$$[n-1]$$$ such that $$$|\sigma(i) - \sigma(i+1)| = 1$$$ for all $$$i \in S$$$. Then we can see that the number of permutations is determined by the following: First arbitrarily permute order the contiguous blocks of $$$[n]$$$ determined by $$$S$$$, then for each block of length $$$\geq 2$$$, it can be ascending or descending, this contributes a factor of $$$2$$$. Therefore, we could write down the generating function of the sequence $$$a_n$$$ as

$$$ \begin{align*} \sum_{n\geq 0} a_n x^n &= \left(\sum_{n\geq 0} n! T^n\right) \circ \left(x - 2 x^2 + 2x^3 - 2x^4 + \cdots\right)\\ &= \left(\sum_{n\geq 0} n! T^n\right) \circ \left(x\frac{1-x}{1+x}\right). \end{align*} $$$

Let $$$A(x)$$$ be the generating function of $$$a_n$$$, if we want to derive a P-recursive relation for $$$a_n$$$, this is equivalent to finding the D-Finite relation for $$$A(x)$$$. The idea is to first analyze the generating function of $$$n! T^n$$$: let

$$$ S(T) = \sum_{n\geq 0} n!T^n, $$$

let $$$s_n = n!$$$ be the sequence of $$$S(T)$$$, then by the recurrence relation

$$$ s_n = n s_{n-1} + [n=0], $$$

we can turn this into a D-finite relation

$$$ S(T) = (S\cdot T)' \cdot T +1 = 1 + TS + T^2S'. $$$

Therefore, since $$$A(x) = S(g(x))$$$, where $$$g(x) = x\frac{1-x}{1+x}$$$, we can derive the D-Finite relation for $$$A(x)$$$:

$$$ A(x) = S \circ g = 1 + gA + g^2 S'\circ g. $$$

recall that $$$A' = g' S'(g)$$$, we can plug in $$$S'\circ g = A' / g'$$$, and we can get the D-Finite relation for $$$A(x)$$$:

$$$ A = 1 + gA + \frac{g^2 A'}{g'}. $$$

After some simplification, we can get the desired P-recursive relation for $$$a_n$$$. I omit the computation details here.

The derivation above is a typical example of the application of the following basic property of D-Finite functions:

Let $$$f$$$ be a D-Finite function, and $$$q$$$ be a rational function, then $$$f\circ q$$$ is D-Finite.

This leads to a fundamental reason why D-Finite functions are useful: They satisfy so many nice closure properties, and these properties are almost complete for explaining the ubiquity of P-recursive sequences in combinatorics.

Multivariate D-Finite Functions

We are already more-or-less familiar with one variable D-Finite functions, whose associated sequence is P-recursive (or called holonomic). A natural question is to ask the correct generalization of D-Finite functions to multiple variables.

The correct definition is due to Richard P. Stanley.

Def 1. Let $$$K$$$ be a field of characteristic zero. A formal power series $$$f \in K[ [ X_1,\dots,X_d ] ]$$$ is called D-finite if the $$$K(X_1,\dots,X_d)$$$-vector space spanned by all derivatives of $$$f$$$ is finite-dimensional.

A caveat is that for multiple variables, the description of the property of the associated sequence of a D-Finite series, is not as nice as in the univariate case. However, we will see why the definition is still useful later.

Similar to the univariate case, multivariate D-Finite functions have nice closure properties.

Theorem. If $$$f, g \in K[ [ X_1,\dots,X_d ] ]$$$ is D-Finite, then

For any $$$a,b \in K(X_1,\dots,X_d)$$$, if $$$af+bg$$$ is still a formal power series, it is D-Finite.
The product $$$fg$$$ is D-Finite.
If $$$u_1,\dots,u_d \in K[ [ Y_1,\dots,Y_e ] ]$$$ is algebraic, and the composition $$$f(u_1,\dots,u_d)$$$ is well-defined in some sense, then it is D-Finite.
(Lipshitz, 1988) If $$$f$$$ is D-Finite, then the diagonal operator $$$\Delta_{X_{d-1},X_d}\colon K[ [ X_1,\dots,X_d ] ] \to K[ [ X_1,\dots,X_{d-2},U ] ]$$$

$$$ \left(\sum_{i_1,\dots,i_d\geq 0} f_{i_1,\dots,i_d} X_1^{i_1}\cdots X_d^{i_d}\right) \mapsto \sum_{i_1,\dots,i_{d-2},j\geq 0} f_{i_1,\dots,i_{d-2},j,j} X_1^{i_1}\cdots X_{d-2}^{i_{d-2}} U^j. $$$

maps D-Finite functions to D-Finite functions.

The first two properties shows that D-Finite functions in $$$K[ [ X_1,\dots,X_d ] ] \otimes_{K[X_1,\dots,X_d]} K(X_1,\dots,X_d)$$$ form a $$$K(X_1,\dots,X_d)$$$-algebra.

The first three properties are easy to prove and the last one needs more insights (but is also the most useful one).

I'll not present the proofs here; interested readers can refer to Manuel Kauers's recent book D-Finite Functions. However, I want to remark that all these properties can be made constructive.

Make Practical

One could prove that Def 1 is equivalent to the following definition:

Def 2. A formal power series $$$f \in K[ [ X_1,\dots,X_d ] ]$$$ is D-Finite if for any $$$1\leq \ell \leq d$$$, there exists $$$k_\ell$$$ such that $$$\partial_{X_\ell}^{k_\ell} f$$$ can be $$$K(X_1,\dots,X_d)$$$-linearly expressed by $$$\partial_{X_\ell}^{j} f$$$ for $$$0\leq j < k_\ell$$$.

Therefore, the constructiveness means that:

We could use data of the linear relation in Def 2 to represent a D-Finite function. (Rigorously speaking, the solution of this system of equations didn't fully determine the D-Finite function, but it is enough for most purposes.)
The above closure properties can be proved by an algorithm that takes the linear relations as input, and outputs the linear relations of the output D-Finite function.

Now we take the simplest example of the sum of two univariate D-Finite functions. Let $$$f,g$$$ be D-Finite, and suppose they satisfy the linear relations

$$$ \partial^n_X f = \sum_{j < n} p_j(X) \partial^j_X f, \quad \partial^m_X g = \sum_{j < m} q_j(X) \partial^j_X g. $$$

Then we can derive the linear relation for $$$f+g$$$: Recursively compute the representation of $$$\partial_X^i (f+g)$$$, in terms of $$$\partial_X^j f$$$ and $$$\partial_X^k g$$$ for $$$j < n$$$ and $$$k < m$$$. Since they form an $$$(n+m)$$$-dimensional vector space, we can always find a linear relation for $$$\partial_X^i (f+g)$$$ between $$$0\leq i\leq n+m$$$.

For example, in jqdai0815's blog, he showed that the power series involved in the problem Chinese Elephant Chess is D-Finite. He suspected that it would be non-practical to compute the sequence via the P-recursive relation, but in fact, I implemented the above algorithm to automatically compute the relation, and then compute the sequence in linear time. Turns out that this code is, in fact, quite efficient. The code is attached below as a reference.

Code

#include <bits/stdc++.h>

#define LOG(FMT...) fprintf(stderr, FMT)

using namespace std;

typedef long long ll;
typedef unsigned long long ull;

// mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());

const uint Mod = 998244353;

constexpr uint norm(uint v) { return v >= Mod ? v - Mod : v; }

struct Z {
  uint v;

  constexpr Z(uint v = 0) : v(v) {}

  constexpr Z(int v) : v(norm(static_cast<uint>(v % Mod + Mod))) {}

  inline constexpr Z operator+(const Z &rhs) const { return norm(v + rhs.v); }

  inline constexpr Z operator-(const Z &rhs) const { return norm(v + Mod - rhs.v); }

  inline constexpr Z operator-() const { return norm(Mod - v); }

  inline constexpr Z operator*(const Z &rhs) const { return static_cast<uint>(static_cast<ull>(v) * rhs.v % Mod); }

  inline constexpr Z inv() const;

  inline constexpr Z operator/(const Z &rhs) const { return *this * rhs.inv(); }

  constexpr operator uint() const { return v; }
};

inline constexpr Z &operator+=(Z &lhs, const Z &rhs) { return lhs = lhs + rhs; }

inline constexpr Z &operator-=(Z &lhs, const Z &rhs) { return lhs = lhs - rhs; }

inline constexpr Z &operator*=(Z &lhs, const Z &rhs) { return lhs = lhs * rhs; }

inline constexpr Z &operator/=(Z &lhs, const Z &rhs) { return lhs = lhs / rhs; }

inline constexpr Z pow(Z base, uint exp) {
  Z res = 1;
  for (; exp; base *= base, exp >>= 1)
    if (exp & 1)
      res *= base;
  return res;
}

inline constexpr Z Z::inv() const { return pow(*this, Mod - 2); }

int nPeak;

struct P : vector<Z> {
  P() : vector(1) {}

  P(unsigned long n) : vector(n) {}

  P(const initializer_list<value_type> &il) : vector(il) {}

  ~P() { nPeak = max(nPeak, (int) size()); }

  int deg() const { return (int) size() - 1; }

  bool trim() {
    int k = size();
    while (k && !operator[](k - 1)) --k;
    resize(k);
    return k;
  }

  P operator-() const {
    P ret(size());
    for (int i = 0; i < size(); ++i) ret[i] = -operator[](i);
    return ret;
  }
};

P operator*(const P &a, const P &b) {
  int n = a.deg(), m = b.deg();
  if (n == -1 || m == -1) return P();
  P c(n + m + 1);
  for (int i = 0; i <= n + m; ++i)
    for (int j = max(0, i - m); j <= min(i, n); ++j)
      c[i] += a[j] * b[i - j];
  return c;
}

P operator*(const P &a, const Z &z) {
  P c(a);
  for (Z &x : c) x *= z;
  return c;
}

P operator+(const P &a, const P &b) {
  if (a.size() >= b.size()) {
    P c = a;
    for (int i = 0; i < b.size(); ++i) c[i] += b[i];
    return c;
  }
  P c = b;
  for (int i = 0; i < a.size(); ++i) c[i] += a[i];
  return c;
}

P operator-(const P &a, const P &b) { return a + -b; }

bool operator==(const P &a, const P &b) {
  if (a.size() < b.size())
    return equal(a.begin(), a.end(), b.begin()) && count(b.begin() + a.size(), b.end(), Z(0)) == b.size() - a.size();
  return equal(b.begin(), b.end(), a.begin()) && count(a.begin() + b.size(), a.end(), Z(0)) == a.size() - b.size();
}

P gcd(P a, P b) {
  if (!a.trim()) return b;
  if (!b.trim()) return a;
  if (a.size() < b.size()) swap(a, b);
  while (b.trim()) {
    Z in = b.back().inv();
    for (Z &x : b) x *= in;
    int n = a.deg(), m = b.deg();
    for (int i = n; i >= m; --i) {
      for (int j = 1; j <= m; ++j) a[i - j] -= a[i] * b[m - j];
      a[i] = 0;
    }
    swap(a, b);
  }
  return a;
}

P div(P a, P b) {
  Z in = b.back().inv();
  for (Z &x : b) x *= in;
  int n = a.deg(), m = b.deg();
  P ret(n - m + 1);
  for (int i = n; i >= m; --i) {
    ret[i - m] = a[i] * b[m];
    for (int j = 1; j <= m; ++j) a[i - j] -= a[i] * b[m - j];
  }
  for (Z &x : ret) x = x * in;
  return ret;
}

struct Q {
  P x, y;

  Q(const P &x = P(), const P &y = {Z(1)}) : x(x), y(y) {}

  Q operator+(const Q &rhs) const { return Q(x * rhs.y + y * rhs.x, y * rhs.y); }

  Q operator-() const { return Q(-x, y); }

  Q operator-(const Q &rhs) const { return *this + -rhs; }

  Q operator*(const Q &rhs) const { return Q(x * rhs.x, y * rhs.y); }

  Q operator*(const Z &rhs) const { return Q(x * rhs, y * rhs); }

  Q inv() const { return Q(y, x); }

  Q operator/(const Q &rhs) const { return *this * rhs.inv(); }

  bool operator==(const Q &rhs) const { return x * rhs.y == y * rhs.x; }

  bool operator!=(const Q &rhs) const { return !operator==(rhs); }

  void simplify() {
    y.trim();
    if (!x.trim()) {
      y = P{Z(1)};
      return;
    }
    P g = gcd(x, y);
    x = div(x, g);
    y = div(y, g);
  }
};

P der(P a) {
  if (a.empty()) return a;
  for (int i = 1; i < a.size(); ++i) a[i] *= i;
  a.erase(a.begin());
  return a;
}

Z eval(const P &a, const Z &z) {
  Z v = 0;
  for (int i = a.deg(); i >= 0; --i)
    v = v * z + a[i];
  return v;
}

vector<Z> invs(const vector<Z>& vec) {
  vector<Z> pre(vec.size()), ret(vec.size());
  pre[0] = 1;
  for (int i = 1; i < vec.size(); ++i)
    pre[i] = pre[i - 1] * vec[i - 1];
  Z tot = accumulate(vec.begin(), vec.end(), Z(1), multiplies<Z>()).inv();
  for (int i = (int)vec.size() - 1; i >= 0; --i) {
    ret[i] = tot * pre[i];
    tot *= vec[i];
  }
  return ret;
}

Q der(const Q &a) {
  return Q(der(a.x) * a.y - der(a.y) * a.x, a.y * a.y);
}

struct Q_Basis {
  int dim, id;
  vector<vector<Q>> basis, augment;

  Q_Basis(int dim) : dim(dim), id(), basis(dim), augment(dim) {}

  vector<Q> insert(vector<Q> vec) {
    vector<Q> rep(dim + 1);
    rep[id++] = Q({Z(1)});
    for (int i = 0; i < dim; ++i) {
      if (vec[i] != Q()) {
        if (basis[i].empty()) {
          for (int j = i + 1; j < dim; ++j) {
            vec[j] = vec[j] / vec[i];
            vec[j].simplify();
          }
          for (int j = 0; j < id; ++j) {
            rep[j] = rep[j] / vec[i];
            rep[j].simplify();
          }
          vec[i] = Q({Z(1)});
          basis[i] = vec;
          augment[i] = rep;
          return {};
        } else {
          for (int j = i + 1; j < dim; ++j)
            vec[j] = vec[j] - vec[i] * basis[i][j];
          for (int j = 0; j < id; ++j)
            rep[j] = rep[j] - vec[i] * augment[i][j];
          vec[i] = Q();
        }
      }
    }
    return rep;
  }
};

// m = size(), sum_m A_m(x)P^(m) = 0
using ODE = vector<P>;
using PRec = vector<P>;

PRec DFinite_genPRec(const ODE &ode) {
  int offset = numeric_limits<int>::max();
  int n = ode.size() - 1, m = numeric_limits<int>::min();
  for (int i = 0; i <= n; ++i) {
    if (ode[i].empty()) continue;
    m = max(m, (int) ode[i].size() - 1 - i);
    int j = 0;
    while (ode[i][j] == 0) ++j;
    offset = min(offset, j - i);
  }
  m -= offset;
  PRec rec(m + 1, P(n + 1));
  P fall{Z(1)};
  for (int i = 0; i <= n; ++i) {
    P coe = fall;
    for (int j = 0; j < (int) ode[i].size() - i - offset; ++j) {
      if (j + i + offset >= 0) rec[j] = rec[j] + coe * ode[i][j + i + offset];
      coe = div(coe * P{-Z(i + j), Z(1)}, P{-Z(j), Z(1)});
    }
    fall = fall * P{-Z(i), Z(1)};
  }
  for (int i = 0; i <= m; ++i) rec[i].trim();
  return rec;
}

struct Eval {
  PRec prec;
  vector<Z> coeff;

  Eval(const PRec &prec) : prec(prec) {}

  int pre(int n) {
    coeff.resize(n + 1);
    for (int i = 0; i <= n; ++i) coeff[i] = eval(prec[0], i);
    int r = 0;
    for (int i = n; i; --i)
      if (coeff[i] == 0) {
        r = i;
        break;
      }
    return r;
  }

  vector<Z> post(vector<Z> init) {
    int start = init.size();
    auto nvs = invs(vector<Z>(coeff.begin() + start, coeff.end()));
    init.resize(coeff.size());
    for (int i = start; i < coeff.size(); ++i) {
      for (int j = 1; j < min(i + 1, (int)prec.size()); ++j)
        init[i] += init[i - j] * eval(prec[j], i);
      init[i] = init[i] * -nvs[i - start];
    }
    return init;
  }
};

ostream &operator<<(ostream &os, const P &p) {
  if (p.empty()) return os << 0;
  os << p[0];
  for (int i = 1; i < p.size(); ++i)
    os << " + " << p[i] << "x^" << i;
  return os;
}

ostream &operator<<(ostream &os, const Q &q) {
  return os << '(' << q.x << ") / (" << q.y << ')';
}

ostream &operator<<(ostream &os, const vector<P> &ode) {
  for (int i = 0; i < ode.size(); ++i) {
    if (i) os << " + ";
    os << '(' << ode[i] << ")P^(" << i << ")";
  }
  return os << " = 0";
}

template<class T>
istream &operator>>(istream &is, vector<T> &v) {
  for (T &x : v)
    is >> x;
  return is;
}

template<class T>
ostream &operator<<(ostream &os, const vector<T> &v) {
  if (!v.empty()) {
    os << v.front();
    for (int i = 1; i < v.size(); ++i)
      os << ' ' << v[i];
  }
  return os;
}

void simplify(ODE &ode) {
  P g = P();
  for (P &x : ode) {
    x.trim();
    g = gcd(g, x);
  }
  for (P &x : ode) if (!x.empty()) x = div(x, g);
}

ODE DFinite_sum(const ODE &op, const ODE &oq) {
  int n = op.size() - 1, m = oq.size() - 1;
  Q_Basis basis(n + m);
  vector<Q> pd(n + 1), qd(m + 1);
  pd[0] = qd[0] = Q({Z(1)});
  for (int dim = 0; dim <= n + m; ++dim) {
    cerr << "INSERT " << pd << " | " << qd << '\n';
    {
      vector<Q> vec(n + m);
      copy(pd.begin(), pd.begin() + n, vec.begin());
      copy(qd.begin(), qd.begin() + m, vec.begin() + n);
      auto ret = basis.insert(vec);
      if (!ret.empty()) {
        cerr << "OK dim = " << dim << "\n";
        ODE ode(dim + 1);
        P prod = {Z(1)};
        for (int i = 0; i < dim; ++i) prod = prod * ret[i].y;
        ode[dim] = prod;
        for (int i = 0; i < dim; ++i) ode[i] = ret[i].x * div(prod, ret[i].y);
        simplify(ode);
        return ode;
      }
    }
    pd[n] = qd[m] = Q();
    for (int j = n - 1; j >= 0; --j) {
      pd[j + 1] = pd[j + 1] + pd[j];
      pd[j] = der(pd[j]);
    }
    for (int j = m - 1; j >= 0; --j) {
      qd[j + 1] = qd[j + 1] + qd[j];
      qd[j] = der(qd[j]);
    }
    for (int j = 0; j < n; ++j) {
      pd[j] = pd[j] - pd[n] * op[j] / op[n];
      pd[j].simplify();
    }
    for (int j = 0; j < m; ++j) {
      qd[j] = qd[j] - qd[m] * oq[j] / oq[m];
      qd[j].simplify();
    }
  }
  cerr << "!FAILED\n";
  assert(false);
}

ODE DFinite_prod(const ODE &op, const ODE &oq) {
  int n = op.size() - 1, m = oq.size() - 1;
  Q_Basis basis(n * m);
  vector<vector<Q>> p(n + 1, vector<Q>(m + 1));
  p[0][0] = Q({Z(1)});
  for (int dim = 0; dim <= n * m; ++dim) {
    {
      vector<Q> vec(n * m);
      for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j)
          vec[i * m + j] = p[i][j];
      cerr << "INSERT " << vec << '\n';
      auto ret = basis.insert(vec);
      if (!ret.empty()) {
        cerr << "OK dim = " << dim << "\n";
        ODE ode(dim + 1);
        P prod = {Z(1)};
        for (int i = 0; i < dim; ++i) prod = prod * ret[i].y;
        ode[dim] = prod;
        for (int i = 0; i < dim; ++i) ode[i] = ret[i].x * div(prod, ret[i].y);
        simplify(ode);
        return ode;
      }
    }
    for (int i = 0; i < n; ++i) p[i][m] = Q();
    for (int j = 0; j < m; ++j) p[n][j] = Q();
    for (int i = n - 1; i >= 0; --i)
      for (int j = m - 1; j >= 0; --j) {
        p[i + 1][j] = p[i + 1][j] + p[i][j];
        p[i][j + 1] = p[i][j + 1] + p[i][j];
        p[i][j] = der(p[i][j]);
      }
    for (int i = 0; i < n; ++i)
      for (int j = 0; j < m; ++j) {
        p[i][j] = p[i][j] - p[n][j] * op[i] / op[n] - p[i][m] * oq[j] / oq[m];
        p[i][j].simplify();
      }
  }
  cerr << "!FAILED\n";
  assert(false);
}

ODE DFinite_compQ(const ODE &op, Q q) {
  q.simplify();
  int n = op.size() - 1;
  vector<vector<Q>> tri(n + 1, vector<Q>(n + 1));
  tri[0][0] = Q({Z(1)});
  Q d = der(q);
  d.simplify();
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j <= i; ++j) {
      tri[i + 1][j] = tri[i + 1][j] + der(tri[i][j]);
      tri[i + 1][j + 1] = tri[i][j] * d;
    }
    for (int j = 0; j <= i + 1; ++j) tri[i + 1][j].simplify();
  }
  vector<Q> vec(n + 1);
  for (int i = 0; i <= n; ++i) {
    for (int j = op[i].deg(); j >= 0; --j)
      vec[i] = vec[i] * q + Q({op[i][j]});
    vec[i].simplify();
  }
  for (int i = n; i >= 0; --i) {
    cerr << tri[i][i] << '\n';
    vec[i] = vec[i] / tri[i][i];
    vec[i].simplify();
    for (int j = 0; j < i; ++j) {
      vec[j] = vec[j] - vec[i] * tri[i][j];
    }
  }
  P prod = P{Z(1)};
  for (int i = 0; i <= n; ++i) prod = prod * vec[i].y;
  ODE ret(n + 1);
  for (int i = 0; i <= n; ++i) ret[i] = vec[i].x * div(prod, vec[i].y);
  simplify(ret);
  return ret;
}

vector<Z> apply(vector<Z> vec, const ODE &ode) {
  vector<Z> ret(vec.size());
  for (int i = 0; i < ode.size(); ++i) {
    for (int j = 0; j < ode[i].size(); ++j)
      for (int k = j; k < vec.size(); ++k)
        ret[k] += vec[k - j] * ode[i][j];
    for (int j = 1; j < vec.size(); ++j)
      vec[j - 1] = vec[j] * Z(j);
    vec.back() = 0;
  }
  return ret;
}

const ODE ODE_EXP = {{-Z(1)},
                     {Z(1)}};
const ODE ODE_LN = {{Z()},
                    {-Z(1)},
                    {Z(1), -Z(1)}};
const ODE ODE_FACT = {{Z(1)},
                      {-Z(1), Z(3)},
                      {Z(),   Z(), Z(1)}};

const Q A1x({Z(), Z(1), Z(1)}, {Z(2), -Z(2)});
const Q B1x({Z(2), -Z(1)}, {Z(2), -Z(2)});

ODE ode_power(const Z &k) {
  return ODE{{-Z(k)},
             {Z(), Z(1)}};
}

ODE ode_bessel(const Z &k) {
  return ODE({{-Z(1)},
              {k + Z(1)},
              {Z(), Z(1)}});
}

int main() {
#ifdef ELEGIA
  freopen("test.in", "r", stdin);
  int nol_cl = clock();
#endif
  ios::sync_with_stdio(false);
  cin.tie(nullptr);

  int n, m;
  cin >> n >> m;
  if (n > m) swap(n, m);
  Z k = m - n;
  ODE A = DFinite_prod(DFinite_compQ(ODE_EXP, A1x), DFinite_compQ(ode_power(-Z(2).inv()), Q({Z(1), -Z(1)})));
  cerr << "A: " << A << '\n';
  ODE B = DFinite_compQ(ode_bessel(k), B1x * B1x * Q({Z(), Z(1)}));
  cerr << "B: " << B << '\n';
  ODE C = DFinite_compQ(ode_power(k), B1x);
  cerr << "C: " << C << '\n';

  ODE tot = DFinite_prod(B, DFinite_prod(A, C));
  cerr << "TOT ODE: " << tot << '\n';

  PRec prec = DFinite_genPRec(tot);
  cerr << prec << '\n';
  Eval eval(prec);
  cerr << eval.coeff << '\n';
  eval.pre(n);
  auto res = eval.post({Z(1)});
  Z ans = res[n];
  Z v = 1;
  for (int i = 1; i <= k; ++i) v *= i; ans /= v;
  for (int i = 1; i <= n; ++i) ans *= i;
  for (int i = 1; i <= n + k; ++i) ans *= i;
  cout << ans << '\n';

  cerr << "peak size = " << nPeak << '\n';

#ifdef ELEGIA
  LOG("Time: %dms\n", int ((clock()
          -nol_cl) / (double)CLOCKS_PER_SEC * 1000));
#endif
  return 0;
}

More Implications

Then, let me demonstrate a more interesting example, this might explain why D-Finite functions are everywhere.

We are often concerned with a summation involving binomials, and we want to simplify the expression. For example, the following sum

$$$ a_n = \sum_k (-1)^k \binom{n}{k}^3. $$$

Well, this actually can be simplified and is known as Dixon's identity, which has a lot of beautiful proofs.

But in a computational perspective, suppose I just want to compute the sequence $$$a_n$$$ efficiently,

do we need to know such a beautiful identity?

The answer is no. It's very general that such kind of sequences are P-recursive, here is a proof:

We first write down the generating function of the binomial:

$$$ Q(X,Y)= \sum_{n,k} \binom{n}{k} X^n Y^k = \sum_n X^n (1+Y)^n = \frac 1{1-X-XY}, $$$

which is clearly D-Finite. Then consider the generating function

$$$ Q(X_1,Y_1)Q(X_2,Y_2)Q(X_3,Y_3), $$$

is clearly D-Finite. Then we use the diagonal operator to glue the variables $$$X_1,X_2,X_3$$$ together, and $$$Y_1,Y_2,Y_3$$$ together, and we get a D-Finite function

$$$ R(X,Y) = \sum_{n,k}\binom{n}{k}^3 X^n Y^k. $$$

Finally, we plug in $$$Y=-1$$$, and we get the generating function of $$$a_n$$$:

$$$ R(X,-1) = \sum_n a_n X^n. $$$

Since $$$R(X,-1)$$$ is D-Finite, $$$a_n$$$ is P-recursive.

Therefore, we could compute $$$a_n$$$ in linear time by just using the recurrence, or even in $$$O(\sqrt n \log n)$$$ time by the fast evaluation algorithm of P-recursive sequences. Knowing the simplification of the expression doesn't help reducing the time complexity in this sense.

This argument has a vast generalization. In the paper Multiple Binomial Sums, the authors formalized a class of summation involving binomials that captures almost all the known identities. One can mimic the above argument to prove that all these sequences are D-Finite. (In the paper they proved some stronger characterization.)

Then here comes another question:

If we just want to solve a problem in competition, do we need to know the proof of P-Recursiveness?

The answer is also no! If a priori we know that the sequence is P-recursive, we can just compute the first several terms in some brute-force way, and then use Gaussian elimination to find the P-recursive relation. This is usually called Min25-BM since the BM algorithm tells that linear recurrence sequences can be reconstructed efficiently.

Remark: The problem of finding the P-recursive relation can be reduced to the Hermite-Padé approximation problem. While Gaussian elimination takes $$$O(n^3)$$$ time, the Hermite-Padé approximation can be solved $$$O(n^2)$$$ time or even $$$O(n\log^2 n)$$$ time. (I'm a little bit hand-waving here about what $$$n$$$ actually means), though this might not be very useful if one just wants to find a P-recursive relation of constant size.

Elegia's blog

Appetizer

Multivariate D-Finite Functions

Make Practical

More Implications