# Faster multivariate integration in D-modules

Hadrien Brochet, Frédéric Chyzak, Pierre Lairez

► To cite this version:

Hadrien Brochet, Frédéric Chyzak, Pierre Lairez. Faster multivariate integration in D-modules. 2025. [https://hal.science/hal-05289785v1](https://hal.science/hal-05289785v1)

HAL Id: hal-05289785

[https://hal.science/hal-05289785v1](https://hal.science/hal-05289785v1)

Preprint submitted on 30 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Faster multivariate integration in D-modules

Hadrien Brochet¹, Frédéric Chyzak¹, and Pierre Lairez¹

¹Inria, 91120, Palaiseau, France

April 17, 2025

## Abstract

We present a new algorithm for solving the reduction problem in the context of holonomic integrals, which in turn provides an approach to integration with parameters. Our method extends the Griffiths-Dwork reduction technique to holonomic systems and is implemented in Julia. While not yet outperforming creative telescoping in D-finite cases, it enhances computational capabilities within the holonomic framework. As an application, we derive a previously unattainable differential equation for the generating series of 8-regular graphs.

## 1 Introduction

Symbolic integration is a fundamental problem in computer algebra, with deep connections to combinatorics, number theory, and mathematical physics. In the vast landscape of integrable functions, holonomic functions—those multivariate functions satisfying sufficiently many independent linear partial differential equations with polynomial coefficients—form a particularly rich and structured class. This class includes many special functions of mathematical and physical interest, such as exponential functions, logarithms, polylogarithms, elliptic integrals and various hypergeometric functions. In the realm of univariate functions, the classical problem of integrating elementary functions in terms of elementary functions does not always have a solution—as governed by Liouville’s theorem. Holonomic functions, in contrast, exhibit a different structure that remains closed under integration, including under definite integration depending on parameters. In the holonomic setting (and its variants), symbolic integration revolves around two problems: integration with a parameter and reduction.

**Integration with a parameter.** Given a function $f(t, x_1, \dots, x_n)$ of $1+n$ variables satisfying a suitable system of linear partial differential equations with polynomial coefficients, we aim at computing a linear differential equation satisfied by the integral

$$I(t) = \int_D f(t, x_1, \dots, x_n)dx_1 \cdots dx_n.$$
To obtain an algebraic formulation of the problem and eliminate the analytic aspects, we need some simplifying assumptions. We consider a function space M in which f lies such that:

(a) M is closed under differentiation by t and the $x_i$, and by multiplication by t and the $x_i$;

(b) $g \mapsto \int_D g(x_1, \dots, x_n)dx_1 \dots dx_n$ is well defined on M and commutes with differentiation with respect to t and multiplication by t;

(c) For any $g \in M$ and any $1 \le i \le n$, $\int_D \frac{\partial}{\partial x_i}g(x_1, \dots, x_n)dx_1 \dots dx_n = 0$.

Then the problem reduces to finding polynomials $p_0(t), \dots, p_r(t)$, with $p_r$ nonzero, and functions $g_1, \dots, g_n \in M$ such that

$$p_0(t)f + p_1(t)\frac{\partial f}{\partial t} + \dots + p_r(t)\frac{\partial^r f}{\partial t^r} = \frac{\partial g_1}{\partial x_1} + \dots + \frac{\partial g_n}{\partial x_n} \quad (1)$$

Indeed, after integrating both sides over D, the hypotheses above imply that

$$p_0(t)I + p_1(t)I' + \dots + p_r(t)I^{(r)} = 0, \quad (2)$$

which is the kind of relations we aim to compute.

To make one step further towards algebra, we introduce the Weyl algebra $W_{t,x}$, which is the non-commutative algebra generated by $t, x_1, \dots, x_n, \partial_t, \partial_{x_1}, \dots, \partial_{x_n}$, and the usual relations $uv = vu, \partial_u \partial_v = \partial_v \partial_u$ and $\partial_u u = u \partial_u + 1$, for any distinct $u, v \in \{t, x_1, \dots, x_n\}$. Hypothesis (a) translates to the statement that M is a left module over $W_{t,x}$. To ensure that (1) has a solution, we require more specifically that:

(d) M is a holonomic $W_{t,x}$-module.

The concept of holonomy embodies the idea of a function satisfying “sufficiently many independent linear PDEs”, see Section 2.2 for more details.

From the algorithmic point of view, we want an algorithm that takes as input a description of M as a $W_{t,x}$-module, with generators and relations, and computes polynomials $p_0(t), \dots, p_r(t)$ such that (1) holds for some $g_1, \dots, g_n$ in M, which we usually do not need to compute.

**The reduction problem.** In absence of a parameter, integrals are constants and holonomic methods do not directly compute them. However, holonomic methods can be used to find relations between integrals. Consider the Weyl algebra $W_x$ in the variables $x_1, \dots, x_n$ and a holonomic $W_x$-module M. As above, we algebraically interpret integration as a linear map on M that vanishes on the subspace $\partial_{x_1}M + \dots + \partial_{x_n}M$, which we denote simply $\partial M$ below. So, finding a relation between the integrals of functions $f_1, \dots, f_r \in M$ means finding constants $c_1, \dots, c_r$, with $c_r \ne 0$, such that

$$c_1 f_1 + \dots + c_r f_r \in \partial M.$$
Computing normal forms in the quotient $M/\partial M$ is one way to find this sort of relations. Owing to the assumption of the holonomy of $M$, a classical result is that this quotient is finite-dimensional over the base field [7, Theorem 6.1 of Chapter 1, combined with the example that precedes it, for $p = 2n$]. It is now well understood that the integration problem may be tackled through the reduction problem [9, 10, 17, 18, 16, 11, 33], in a way that departs from algorithms based on earlier approaches, namely by a linear ansatz [52, 36] or by following Zeilberger's method [53, 20].

**Contributions.** In the holonomic context, we propose a new algorithm for the problem of reduction (Section 3). In some aspects, this is a generalization to holonomic systems of the Griffiths–Dwork reduction method for homogeneous rational functions [32, 12, 39] (Section 3.4). This algorithm yields a new algorithm for integration with a parameter (Section 4). We provide a Julia [6] implementation of our algorithm (Section 5). Although the new algorithm is still not on par with best implementations following Zeilberger's approach and generalizations on their home turf (i.e., D-finite functions, see below), it improves the state of the art in the holonomic context. As an application, we were able to compute a differential equation for the generating series of 8-regular graphs, for which Zeilberger's approach is theoretically not suited, and which was also previously unattainable by dedicated methods (Section 6).

**D-finiteness versus holonomy.** The creative telescoping approach to symbolic integration (see [21] for a review) relies on *D-finiteness* instead of holonomy. Instead of working with the Weyl algebra $W_{t,x}$ and with holonomic $W_{t,x}$-modules, this approach considers the rational Weyl algebra $W_{t,x}(t, \mathbf{x})$, which is the Weyl algebra extended with rational functions in $t$ and $\mathbf{x}$, and $W_{t,x}(t, \mathbf{x})$-modules that are finite-dimensional over $\mathbb{K}(t, \mathbf{x})$, where $\mathbb{K}$ is the base field. The nuance has deep concrete implications.

Expressivity is an argument in favor of holonomy: we can express a wider class of integrals with holonomy than with D-finiteness. Integrals over semi-algebraic sets is a prominent example [43]. In general, it is always possible to construct a D-finite module from a holonomic module $M$: it is enough to consider the localization $\mathbb{K}(t, \mathbf{x}) \otimes_{\mathbb{K}[t, \mathbf{x}]} M$, but this operation may lose important information. For example, it is possible that $\mathbb{K}(t, \mathbf{x}) \otimes_{\mathbb{K}[t, \mathbf{x}]} M = 0$, as it happens when enumerating $k$-regular graphs (see Section 6), making it impossible to apply any creative-telescoping algorithm over rational functions in a relevant way.

As for efficiency, approaches based on D-finiteness (implemented in Mathematica [36, 37], Maple [20, 11], and Sagemath [35]) are far superior, when they apply, to that based on holonomy [44] (implemented in Macaulay2, Singular, and Risa/Asir). Understanding and bridging this gap to achieve both efficiency and expressivity remains a significant challenge, and this paper represents a step forward in that direction.

In this paper, we consider a mixed approach by using modules over $W_{t,x}(t)$,
the Weyl algebra extended with rational functions in $t$ only: this amounts to considering holonomy with respect to the $x$ and D-finiteness with respect to the parameter $t$. This enables the use of reductions over the base field $K(t)$ to compute integrals with a parameter.

**Related work.** In the context of D-finiteness, the problem of integration with parameters is addressed by algorithms of Chyzak [20] and Koutschan [36]. A more recent research line addresses the integration problem by solving first the reduction problem [9, 10, 17, 18, 16, 11, 33, 15]. In a holonomic context, the integration problem and the reduction problem have been addressed by Takayama [50] and Oaku and Takayama [44], without making efficiency their main goal. In this work, we forgo the minimality of the order of the output relation (2) in order to bypass the machinery related to $b$-functions. Already Takayama’s algorithm [50] made a similar compromise. However, in comparison to this algorithm, we leverage a Gröbner basis technique to obtain a first reduction, which has a lot of structure and can be computed efficiently, but which is not enough to detect all relations between integrals. This is completed by another reduction, more direct and less structured, which happens in a smaller dimensional space than what would arise with Takayama’s algorithm, thanks to the first reduction.

In the case of integrals of the form

$$ \int a(x_1, \dots, x_n) \exp(f(x_1, \dots, x_n)) dx_0 \dots dx_n, $$

where $a(x_1, \dots, x_n)$ is a polynomial and $f(x_1, \dots, x_n)$ a homogeneous polynomial, the first reduction echoes the Griffiths–Dwork reduction [32, 12], while the second reduction echoes Lairez’s reduction algorithm in [38]. In the similar context of rational integrals, which are of great importance in the computation of Feynman integrals, the reduction is addressed by Laporta [41].

In the context of the combinatorics of $k$-regular graphs, first algorithms for computing linear differential equations satisfied by their counting generating functions were developed by Chyzak, Mishna, and Salvy [23]. This was following works by the combinatorialist Gessel in the 1980s, who introduced a representation of the generating functions as a scalar product in the theory of symmetric functions. A faster method was very recently introduced by Chyzak and Mishna [22], based on the same scalar-product representation but following an approach reminiscent of reductions. This was the starting point of our interest, making us rethink the representation to have the algorithms of the present work apply directly to the problem.

# 2 Computing with Weyl algebras

## 2.1 Weyl algebras

Let $K$ be a field of characteristic zero, typically $\mathbb{Q}$ or $\mathbb{Q}(t)$. Let $W_x$ denote the $n$th Weyl algebra $K[x] \langle \partial_x \rangle$ with generators $x = (x_1, \dots, x_n)$ and $\partial_x = (\partial_1, \dots, \partial_n)$,
and relations $\partial_i x_i = x_i \partial_i + 1$, $x_i x_j = x_j x_i$, $\partial_i \partial_j = \partial_j \partial_i$ and $x_i \partial_j = \partial_j x_i$ whenever $i \neq j$. We refer to [24, Chapters 1–10] for a complete introduction to these algebras covering most needs of the present article, or to [8, Chapter 5] for a denser alternative. We often need to highlight one variable with a specific role, in which case we use the name $t$ for the distinguished variable. Correspondingly, we will write $W_{t,\mathbf{x}}$ for the $(1+n)$th Weyl algebra, and we will write $W_t$ for the special case $n=0$. We also define $W_{t,\mathbf{x}}(t)$ as the algebra $K(t) \otimes_{K[t]} W_{t,\mathbf{x}}$ where the variable $t$ is rational and the variables $\mathbf{x}$ are polynomial. For non-zero $r \in \mathbb{N}$, we also consider Cartesian powers of these algebras, $W_{\mathbf{x}}^r$, $W_{\mathbf{x}}(t)^r$, etc., which we view as modules over $W_{\mathbf{x}}$ or $W_{\mathbf{x}}(t)$, as relevant. Each element of the module $W_{\mathbf{x}}^r$ decomposes uniquely in the basis of $K$-vector space

$$M_{\mathbf{x},r} = \{ \mathbf{x}^{\alpha} \partial_{\mathbf{x}}^{\beta} e_i \mid \alpha, \beta \in \mathbb{N}^n, i \in \{1, \dots, r\} \}.$$

Given an element $p = \sum_{\alpha, \beta, i} a_{\alpha, \beta, i} \mathbf{x}^{\alpha} \partial^{\beta} e_i$ of $W_{t,\mathbf{x}}^r$, we define the degree of $p$ as

$$\deg(p) = \max\{|\alpha| + |\beta| \mid \exists i \in \{1, \dots, r\}, a_{\alpha, \beta, i} \neq 0\},$$

where $|\alpha|$ and $|\beta|$ denote the sums $\alpha_1 + \dots + \alpha_n$ and $\beta_1 + \dots + \beta_n$. Definitions for the algebra $W_{\mathbf{x}}$ mimic the case $r=1$, just not considering any $e_i$, and definitions for the module $W_{t,\mathbf{x}}^r$ are just a special notation when $n$ is replaced with $n+1$: a vector basis of $W_{\mathbf{x}}$ is

$$M_{\mathbf{x}} = \{ \mathbf{x}^{\alpha} \partial_{\mathbf{x}}^{\beta} \mid \alpha, \beta \in \mathbb{N}^n \},$$

and given an element $p = \sum_{\alpha, \beta} a_{\alpha, \beta} \mathbf{x}^{\alpha} \partial_{\mathbf{x}}^{\beta}$ of $W_{\mathbf{x}}$, its degree is

$$\deg(p) = \max\{|\alpha| + |\beta| \mid a_{\alpha, \beta} \neq 0\};$$

a vector basis of $W_{t,\mathbf{x}}^r$ is

$$M_{t,\mathbf{x},r} = \{ t^{\alpha} \mathbf{x}^{\beta} \partial_{t}^{\gamma} \partial_{\mathbf{x}}^{\delta} e_i \mid \alpha, \gamma \in \mathbb{N}, \beta, \delta \in \mathbb{N}^n, i \in \{1, \dots, r\} \},$$

and given an element $p = \sum_{\alpha, \beta, \gamma, \delta, i} a_{\alpha, \beta, \gamma, \delta, i} t^{\alpha} \mathbf{x}^{\beta} \partial_{t}^{\gamma} \partial_{\mathbf{x}}^{\delta} e_i$ of $W_{t,\mathbf{x}}^r$, its degree is

$$\deg(p) = \max\{\alpha + |\beta| + \gamma + |\delta| \mid \exists i \in \{1, \dots, r\}, a_{\alpha, \beta, \gamma, \delta, i} \neq 0\}.$$

## 2.2 Holonomic modules

Let $S$ be a submodule of $W_{\mathbf{x}}^r$. We recall the classical definition of a holonomic $W_{\mathbf{x}}$-module by means of the dimension of the quotient module $M = W_{\mathbf{x}}^r/S$. We point out that any $W_{\mathbf{x}}$-module of finite type is isomorphic to a module of this form. The *Bernstein filtration* [4] of the algebra $W_{\mathbf{x}}$ is the sequence of $K$-vector spaces $\mathcal{F}_m$ defined by

$$\mathcal{F}_m = \{ P \in W_{\mathbf{x}} \mid \deg(P) \leq m \}.$$

A filtration of the module $M$ that is adapted to $(\mathcal{F}_m)_{m \geq 0}$ is the sequence of $K$-linear subspaces $\Phi_m \subseteq M$ defined by

$$\Phi_m = \text{image in } M \text{ of } \{ P \cdot e_i \mid P \in \mathcal{F}_m, 1 \leq i \leq r \}.$$
Those filtrations are compatible with the algebra and module structures: both $\mathcal{F}_m\mathcal{F}_{m'} \subseteq \mathcal{F}_{m+m'}$ and $\mathcal{F}_m\Phi_{m'} \subseteq \Phi_{m+m'}$ hold. There exists a polynomial $p \in \mathbb{K}[m]$ called the *Hilbert polynomial* of $M$ that satisfies $\dim_{\mathbb{K}}(\Phi_k) = p(k)$ for any sufficiently large $k$. The *dimension* of the module $M$ is the degree $d$ of the polynomial $p$. The integer $d$ clearly lies between 0 and $2n$. It was proved by Bernstein that if $M$ is non-zero, then $d$ is larger than or equal to $n$ [24, Theorem 9.4.2]. When the dimension of $M$ is exactly $n$ or when $M$ is the zero module, we say that the module $M$ is *holonomic*. (Here, we follow the tradition in [8, Chapter 5] and in [24] to consider the zero module as holonomic. By way of comparison, [7] speaks of a module “in the Bernstein class” to refer to a non-zero holonomic module.)

## 2.3 Gröbner bases in Weyl algebras and their modules

Despite their non-commutative nature, by a *monomial* we will mean an element of the vector bases $M_{\mathbf{x},r}$, $M_{\mathbf{x}}$, and $M_{t,\mathbf{x},r}$. A *monomial order* $\prec$ on $W_{\mathbf{x}}^r$ is a well-ordering on $M_{\mathbf{x},r}$ that satisfies for any $i,j \in \{1, \dots, r\}$ and any exponents $\alpha, \beta, \alpha_1, \beta_1, \alpha_2, \beta_2 \in \mathbb{N}^n$

$$ \mathbf{x}^{\alpha_1} \partial_{\mathbf{x}}^{\beta_1} e_i \prec \mathbf{x}^{\alpha_2} \partial_{\mathbf{x}}^{\beta_2} e_j \implies \mathbf{x}^{\alpha_1+\alpha} \partial_{\mathbf{x}}^{\beta_1+\beta} e_i \prec \mathbf{x}^{\alpha_2+\alpha} \partial_{\mathbf{x}}^{\beta_2+\beta} e_j. $$

Given an operator $P \in W_{\mathbf{x}}^r$, we define its *support* $\text{supp}(P)$ as the set of all monomials appearing with non-zero coefficient in the decomposition of $P$ with coefficients on the left in the basis $M_{\mathbf{x},r}$. We then define its *leading monomial* $\text{lm}(P)$ as the largest monomial for $\prec$ in $\text{supp}(P)$, its *leading coefficient* $\text{lc}(P)$ as the coefficient of $\text{lm}(P)$ in this decomposition, and its *leading term* $\text{lt}(P)$ as $\text{lc}(P) \text{lm}(P)$. We stress that our definition of a leading monomial makes $\text{lm}(P)$ an element of $W_{\mathbf{x}}^r$, whereas some authors choose to see leading monomials as commutative objects in an auxiliary commutative polynomial algebra, introducing commutative variables $\xi_i$ to replace the $\partial_i$ in monomials. For example, we have $\text{lm}(\partial_1 x_1 e_i) = x_1 \partial_1 e_i$. Note that an essentially equivalent theory could be developed by choosing monomials as elements of the basis consisting of the products $\partial_{\mathbf{x}}^{\beta} \mathbf{x}^{\alpha} e_i$, instead of monomials in $M_{\mathbf{x},r}$.

Computations in Weyl algebras rely heavily on a non-commutative generalization of Gröbner bases. After the original introduction [14, 13, 30] there have been a number of presentations of such a theory, including [48, 34, 42]. A first textbook presentation is [46, Chapter 1]. A recent and simpler introduction can be found in [2]. We now adapt this to our need. A Gröbner basis of a left (resp. right) ideal $I$ of $W_{\mathbf{x}}$ with respect to a monomial order $\prec$ is a finite set $G$ of generators of $I$ such that for any $a \in I$ there exist $g \in G$ and $q \in W_{\mathbf{x}}$ satisfying $\text{lm}(a) = \text{lm}(qg)$ (resp. $\text{lm}(a) = \text{lm}(gq)$). Note that the non-commutativity of the monomials makes them lack divisibility properties that a commutative variant would: we may have $\text{lm}(qg) \neq \text{lm}(q) \text{lm}(g)$ and $\text{lm}(a)$ may not be a multiple of $\text{lm}(g)$, be it on the left or on the right. However, we always have $\text{lm}(a) = \text{lm}(\text{lm}(q) \text{lm}(g))$. (The variables $\xi_i$ introduced by other authors is to avoid this formula.) A Gröbner basis allows to define and compute for
any $a \in W_\mathbf{x}$ a unique representative of $a + I$ in the quotient $W_\mathbf{x}/I$ by means of a non-commutative generalization of polynomial division. We call this unique representative the *remainder of the division of a by the Gröbner basis* $G$ or more shortly *remainder of a modulo the Gröbner basis* $G$. We denote this remainder $\mathrm{LRem}(a, G)$ when $I$ is a left ideal and $\mathrm{RRem}(a, G)$ when $I$ is a right ideal. The concept of Gröbner bases for ideals of $W_\mathbf{x}$ extends to submodules of $W_\mathbf{x}^r$ in the same way as the notion of Gröbner bases for polynomial ideals generalizes to submodules of a polynomial algebra (see [3, Chapter 10.4], [25, Chapter 5], or [1, Section 3.5]). The noetherianity of Weyl algebras implies that any left or right submodule of $W_\mathbf{x}^r$ admits a (finite) Gröbner basis.

Let $A$ be a subvector space of $W_\mathbf{x}^r$. We define $\partial A$ as the $\mathbb{K}$-vector space $\sum_{i=1}^n \partial_i A$. If $A$ is a right $W_\mathbf{x}$-module, then $\partial A$ is also a right $W_\mathbf{x}$-module.

## 2.4 Integration

The integral of a $W_\mathbf{x}$-module $M \simeq W_\mathbf{x}^r/S$ is the $\mathbb{K}$-vector space

$$ M/\partial M \simeq W_\mathbf{x}^r/(S + \partial W_\mathbf{x}^r). \qquad (3) $$

As already mentioned, it is classical that, if $M$ is holonomic, then the integral (3) of $M$ is a finite-dimensional $\mathbb{K}$-linear space [8, Theorem 6.1 of Chapter 1]. Computing relations modulo $\partial M$ is the main matter of this article: given a family in $M$, we want to find a linear dependency relation on its image in the integral module $M/\partial M$, if any such relation exists.

## 2.5 Data structure for holonomic modules

Algorithmically, we only deal with holonomic $W_\mathbf{x}$-modules. They are finitely presented: for such a module $M$, there exist $W_\mathbf{x}$-linear homomorphisms $a$ and $b$ forming an exact sequence

$$ W_{\mathbf{x}}^{s} \xrightarrow{a} W_{\mathbf{x}}^{r} \xrightarrow{b} M \to 0. $$

Equivalently, this means that $M \simeq W_\mathbf{x}^r/S$, where $S$ is the left submodule generated by the image under $a$ of the canonical basis of $W_\mathbf{x}^s$. The module $S$ consists of $W_\mathbf{x}$-linear combinations of the canonical basis of $W_\mathbf{x}^r$, which we always denote $(e_1, \dots, e_r)$. This gives a concrete data structure for representing holonomic modules. It is well known that holonomic modules are *cyclic*, that is generated by a single element [24, Corollary 10.2.6]. This means that we could in principle always assume that $r = 1$. However, some modules have a more natural description with $r > 1$ and transforming the presentation to achieve $r = 1$ has an algorithmic cost that we are not willing to pay. Therefore, we will not assume $r = 1$.

In order to integrate an infinitely differentiable function $f(\mathbf{x})$, we may consider the $W_\mathbf{x}$-module generated by $f$ under the natural action of $W_\mathbf{x}$ on $C^\infty$ functions. Of course, the holonomic approach to symbolic integration will only work
if this module is holonomic. Instead of $W_x \cdot f$, we can also consider any holonomic module that contains it as a submodule.

For example, to integrate a rational function $A/F \in K(\mathbf{x})$, we can consider the module $K[\mathbf{x}][F^{-1}]$, which is holonomic. However, finding a finite presentation $W_x^r/S \simeq K[\mathbf{x}][F^{-1}]$ is not trivial. There are algorithms [45] to solve this problem, but, in terms of efficiency, it is still a practical issue that we do not address in this work. Fortunately, in many cases we can easily construct a holonomic module for integration. See for example Section 6.

# 3 Reductions

We consider the Weyl algebra $W_x$ over a field $K$ and a finitely presented $W_x$-module $M$ given in the form $W_x^r/S$ for some $r \ge 1$ and some submodule $S$ of $W_x^r$ (see Section 2.2). The main objective of the section is to compute normal forms in $M$ modulo $\partial M$, or, equivalently, normal forms in $W_x^r$ modulo $S + \partial W_x^r$. In other words, we want an algorithm that given some $a \in W_x^r$ computes some $[a] \in W_x^r$, and such that $[a] = [b]$ if and only if $a - b \in S + \partial W_x^r$. This goal in only partially reached with a family of reductions $[\cdot]_\eta$ such that for each pair $(a, b)$, there exists $\eta$ such that $[a]_\eta = [b]_\eta$ if and only if $a - b \in S + \partial W_x^r$. The existence of a monomial $\eta$ is not effective, similarly to the maximal total degree to be considered in Takayama's algorithm [49]. This is a step backwards compared to previous methods [44], but computing weaker normal forms allows for more efficient computational methods. Concretely, we do not rely on the computation of $b$-functions.

The present section is organized as follows. In Section 3.1 we define a reduction procedure $[\cdot]$ that partially reduces elements of $W_x^r$ by $S + \partial W_x^r$, in the sense that the procedure will in general not reduce every element of $S + \partial W_x^r$ to zero. In Section 3.2 we define a filtration $(F_{\prec\eta})_{\eta \in M_{x,r}}$ of the vector space $S + \partial W_x^r$ and we give an algorithm to compute a basis of each vector space $[F_{\prec\eta}]$ of reduced forms. Using this basis we define a new reduction $[\cdot]_\eta$ that enhances the first one. In Section 3.3, we provide, for some infinite families $(a_i)_{i \ge 0}$ in $W_x^r$, an algorithm for computing an $\eta$ such that all the $[a_i]_\eta$ lie in a finite-dimensional subspace. In Section 3.4, we consider the case where $S$ is the annihilator of $e^f$ for some homogeneous multivariate polynomial $f$, and we compare our reduction procedures with variants of the Griffiths-Dwork reduction.

## 3.1 Reduction $[\cdot]$ and irreducible elements

**Reduction rules.** Let $\prec$ be a monomial order on $W_x^r$ and let $G$ be a Gröbner basis of $S$ for this order. We define two binary relations $\to_1$ and $\to_2$ on $W_x^r \times W_x^r$ as follows:

* Given $a \in W_x^r$, $\lambda \in K$, $m \in M_x$, and $g \in G$, we write

$$a \to_1 a - \lambda mg$$

if $\mathrm{lm}(mg)$ is in the support of $a$ but not in the support of $a - \lambda mg$.
* Given $a \in W_{\mathbf{x}}^r$, $\lambda \in \mathbb{K}$, $m \in M_{\mathbf{x},r}$, and $i \in \{1, \dots, n\}$, we write

$$a \to_2 a - \lambda \partial_i m$$

if $\text{lm}(\partial_i m)$ is in the support of $a$ but not in the support of $a - \lambda \partial_i m$.

The relation $\to_1$ corresponds to the reduction by the Gröbner basis $G$ of the left module $S$ and the relation $\to_2$ corresponds to the reduction by the Gröbner basis $\{\partial_i e_j \mid i = 1, \dots, n, j = 1, \dots, r\}$ of the right module $\partial W_{\mathbf{x}}^r$. Next, we define $\to$ as the relation $\to_1 \cup \to_2$. That is, $a \to b$ if either $a \to_1 b$ or $a \to_2 b$. The relation $\to^+$ is the transitive closure of $\to$: $a \to^+ b$ if there exist $s \ge 1$ and a sequence of $s$ reductions

$$a \to c_1 \to \dots \to c_s = b \quad (4)$$

for some $c_1, \dots, c_s \in W_{\mathbf{x}}^r$. The relation $\to^*$ is the reflexive closure of $\to^+$: $a \to^* b$ if either $a \to^+ b$ or $a = b$. In this situation we say that $a$ reduces to $b$.

**Irreducible elements.** We say that an element $b$ is irreducible if there is no $c$ such that $b \to c$ and we say that $b$ is a reduced form of $a$ if $b$ is irreducible and $a \to^* b$.

**Lemma 1.** Let $a, b \in W_{\mathbf{x}}^r$. If $a \to^* b$ then $a - b \in S + \partial W_{\mathbf{x}}^r$.

*Proof.* This follows from the definition of $\to_1$ and $\to_2$ since the terms $mg$ and $\partial_i m$ are in $S$ and $\partial W_{\mathbf{x}}^r$, respectively. $\square$

However, the converse of Lemma 1 is not true in general, even when $b = 0$: there may be nonzero irreducible elements in $S + \partial W_{\mathbf{x}}^r$.

**Lemma 2.** The irreducible elements of $W_{\mathbf{x}}^r$ form a $\mathbb{K}$-vector space.

*Proof.* The set $V$ of all irreducible elements contains 0 and is stable by multiplication by $\mathbb{K}$. Let $a, b \in V$ and assume by contradiction that $a + b$ is not irreducible. Then, there exists a monomial $m \in M_{\mathbf{x},r}$ in the support of $a+b$ that can be reduced by $\to$. Because $a+b=b+a$, we can without loss of generality assume that $m$ is in the support of $a$. This contradicts the irreducibility of $a$. Thus $a+b \in V$. $\square$

The vector space of Lemma 2 can be infinite-dimensional, as we now exemplify.

*Example 3.* Let $S = W_{x_1} \partial_1$ be the left ideal generated by $\partial_1$ in the Weyl algebra in one pair of generators, $(x_1, \partial_1)$. Note that $W_{x_1}/S \simeq \mathbb{K}[x_1]$ as $W_{x_1}$-module. Let $\llcorner$ be the lexicographic order for which $\partial_1 \llcorner x_1$. Then, any element of $\mathbb{K}[x_1]$ (as a subspace of $W_{x_1}$) is irreducible.

Irreducible forms can be computed by alternating left reductions with respect to a Gröbner basis of $S$ (representing the rule $\to_1$) and right reductions with respect to a Gröbner basis of $\partial W_{\mathbf{x}}^r$ (representing the rule $\to_2$). This leads to Algorithm 1. Correctness is clear. The algorithm terminates since the largest reducible monomial in $a$, if any, decreases at each iteration of the loop.
<caption>**Algorithm 1** Computation of a reduced form</caption>

**Input:**
* $a \in W_x^r$
* a Gröbner basis $G$ of $S$

**Output:**
* a reduced form of $a$

1 **while** $a$ is not irreducible
2 &nbsp;&nbsp;&nbsp;&nbsp;$a \leftarrow \text{RRem}(a, \{\partial_i e_j \mid i = 1, \dots, n, j = 1, \dots, r\})$
3 &nbsp;&nbsp;&nbsp;&nbsp;$a \leftarrow \text{LRem}(a, G)$
4 **return** $a$

**Definition 4.** We denote by $[a]$ the reduced form of $a \in W_x^r$ that is computed by Algorithm 1.

**Proposition 5.** The map $[\cdot]$ is $\mathbb{K}$-linear.

*Proof.* The maps RRem and LRem are $\mathbb{K}$-linear by the uniqueness of the remainder of a division by a Gröbner basis. Let $\tau(a)$ denote the number of iterations of the while loop in Algorithm 1 on input $a$. Given $\tau \in \mathbb{N}$, let $V_\tau$ denote the set of all $a$ for which $\tau(a) \le \tau$. The restriction of $[\cdot]$ on $V_\tau$ takes the same values as the composition of $\tau$ copies of RRem and $\tau$ copies of LRem in alternation, in which some of the final copies effectively act by the identity as they input irreducible elements. So the restriction of $[\cdot]$ on $V_\tau$ is $\mathbb{K}$-linear as a composition of linear maps. The result follows because $W_x^r = \bigcup_{\tau \ge 0} V_\tau$. $\square$

## 3.2 Computation of the irreducible elements of $S + \partial W_x^r$

Again, we fix an order $\prec$ and a submodule $S$ of $W_x^r$ by considering a Gröbner basis $G$ of it. Let $E$ be the vector space of all irreducible elements of $S + \partial W_x^r$. This vector space can be infinite-dimensional hence we cannot hope to compute all of it. We therefore define a vector-space filtration $(F_{\prec\eta})_{\eta \in M_{x,r}}$ of $S + \partial W_x^r$ by

$$F_{\prec\eta} = \{s+d \in W_x^r \mid s \in S, d \in \partial W_x^r, \text{ and } \max(\mathrm{lm}(s), \mathrm{lm}(d)) \prec \eta\},$$

and a vector-space filtration of $E$ by $E_{\prec\eta} := F_{\prec\eta} \cap E$. We define $F_{\prec\eta}$ and $E_{\prec\eta}$ similarly, by requiring a strict inequality on the maximum of the leading monomials.

Our goal is to obtain an efficient computation of a $\mathbb{K}$-basis of $E_{\prec\eta}$. Let us give an intuitive description of our algorithm. By general properties of Gröbner bases, a non-zero element reduces to zero using the relation $\to_1$ (resp. $\to_2$) if and only if it belongs to $S$ (resp. $\partial W_x^r$). The difficulty arises when both reduction rules can be applied to reduce a monomial. For example, take an element $s$ in $S$ such that $\mathrm{lm}(s) \in \mathrm{lm}(S) \cap \mathrm{lm}(\partial W_x^r)$ and, assuming is can be reduced so as to cancel its leading monomial by using $\to_2$, perform this reduction, that is,
find $s'$ such that $s \to_2 s'$ and $\mathrm{lm}(s') \prec \mathrm{lm}(s)$. In this case, it is possible that $s'$ is neither in $S$ nor in $\partial W_x^r$, making it a good candidate for an element that does not reduce to 0 by $\to$. The following theorem shows more precisely how generators of $E$ can be obtained.

**Theorem 6.** Let $\eta \in M_{\mathbf{x},r}$.

1. If $\eta \notin \mathrm{lm}(S) \cap \mathrm{lm}(\partial W_x^r)$, then $E_{\prec \eta} = E_{\prec \eta}$.

2. If $\eta \in \mathrm{lm}(S) \cap \mathrm{lm}(\partial W_x^r)$, then $E_{\prec \eta} = E_{\prec \eta} + \mathbb{K}a$, for any reduced form $a$ of $mg - \partial_i w$, where $w \in W_x^r$, $m \in W_x$ and $g \in G$ are any elements such that $\eta = \mathrm{lm}(mg) = \mathrm{lm}(\partial_i w)$ and $\mathrm{lc}(mg) = \mathrm{lc}(\partial_i w)$. Moreover, such $m$ and $g$ exist because $G$ is a Gröbner basis of $S$.

*Proof.* For the first point, we prove by contradiction that for any $a \in E_{\prec \eta}$ and any $s \in S$ and $d \in \partial W_x^r$ satisfying $a = s + d$ and $\mathrm{max}(\mathrm{lm}(s), \mathrm{lm}(d)) \prec \eta$, we have in fact $\mathrm{max}(\mathrm{lm}(s), \mathrm{lm}(d)) \prec \eta$. This will imply the equality $E_{\prec \eta} = E_{\prec \eta}$. Let us assume that the equality $\mathrm{max}(\mathrm{lm}(s), \mathrm{lm}(d)) = \eta$ holds. Therefore, either $\mathrm{lm}(s) = \mathrm{lm}(d) = \eta$, or $\mathrm{lm}(s) \prec \mathrm{lm}(d) = \eta$, or $\mathrm{lm}(d) \prec \mathrm{lm}(s) = \eta$. The first case is excluded because we assumed $\eta \notin \mathrm{lm}(S) \cap \mathrm{lm}(\partial W_x^r)$. In both remaining cases it is possible to reduce $\mathrm{lm}(a) = \eta$ with one of the two reduction rules. This contradicts the fact that $a$ is irreducible.

For the second point, let $m, g, w, i$, and $a$ be given as in the statement. We first check that $E_{\prec \eta} + \mathbb{K}a \subseteq E_{\prec \eta}$. It is enough to prove that $a \in E_{\prec \eta}$. By definition, $mg - \partial_i w \in F_{\prec \eta}$, and we check easily that $F_{\prec \eta}$ is stable under $\to$. So $a \in F_{\prec \eta}$. Since $a$ is also irreducible, we have $a \in E_{\prec \eta}$.

Let us prove the other inclusion. Let $f \in E_{\prec \eta}$. Then $f$ is irreducible and of the form $s+d$ for $s \in S$ and $d \in \partial W_x^r$ satisfying $\mathrm{max}(\mathrm{lm}(s), \mathrm{lm}(d)) \prec \eta$. If this inequality is strict, then $f \in E_{\prec \eta}$, proving $f \in E_{\prec \eta} + \mathbb{K}a$. Otherwise, we have the equality $\mathrm{max}(\mathrm{lm}(s), \mathrm{lm}(d)) = \eta$. Let us remark the equality $\mathrm{lm}(s) = \mathrm{lm}(d)$, for otherwise either $\mathrm{lm}(s) \succ \mathrm{lm}(d)$ and $f$ could be reduced using $\to_1$, or $\mathrm{lm}(s) \prec \mathrm{lm}(d)$ and $f$ could be reduced using $\to_2$. So $\eta = \mathrm{lm}(s) = \mathrm{lm}(d)$. This monomial cannot be $\mathrm{lm}(f)$, for otherwise $f$ could be reduced using any of $\to_1$ and $\to_2$. Hence $\mathrm{lm}(f) \prec \eta$ and $\mathrm{lt}(s) = -\mathrm{lt}(d)$. We decompose $s$ and $d$ as $s = \lambda mg + s'$ and $d = -\lambda \partial_i w + d'$ with $\lambda \in \mathbb{K}$, $s' \in S$, $d' \in \partial W_x^r$, and $\mathrm{max}(\mathrm{lm}(s'), \mathrm{lm}(d')) \prec \eta$. Let $h$ denote $mg - \partial_i w$, which, by hypothesis, has $a$ as a reduced form. This implies an equality of the form $h = a + s'' + d''$ with $s'' \in S$, $d'' \in \partial W_x^r$, and $\mathrm{max}(\mathrm{lm}(s''), \mathrm{lm}(d'')) \prec \eta$. We obtain $f = s+d = \lambda(mg - \partial_i w) + s' + d' = \lambda a + b$ with $b = s' + \lambda s'' + d' + \lambda d''$. Since both $f$ and $a$ are irreducible so is $b$, thus $b \in E_{\prec \eta}$, proving that $f$ is in $E_{\prec \eta} + \mathbb{K}a$. $\square$

The meaning of Theorem 6 is that the dimension of the filtration $(E_{\prec \eta})_{\eta}$ is susceptible to increase at $\eta$ only if $\eta \in \mathrm{lm}(S) \cap \mathrm{lm}(\partial W_x^r)$. But this is not necessary as the element $a$ may well be in $E_{\prec \eta}$. The following lemma describes a sufficient condition for this situation.

**Lemma 7.** Let $\eta \in \mathrm{lm}(S) \cap \mathrm{lm}(\partial W_x^r)$. If there exist $g \in G$, $m \in M_x$, and some $i$ such that $\eta = \mathrm{lm}(\partial_i mg)$, then

$$E_{\prec \eta} = E_{\prec \eta}.$$
*Proof.* By Theorem 6, the result reduces to proving that $E_{<\eta}$ contains the reduced form $a$ of some $h = \partial_i mg - \partial_j w$ with $\mathrm{lm}(h) < \eta$. We choose $j = i$ and $w = mg$, so that $h = 0$, which already is irreducible and in $E_{<\eta}$. $\square$

**Corollary 8.** Let $\eta \in M_{x,r}$. Let $H$ be the set of monomials $m \le \eta$ such that $m \in \mathrm{lm}(S) \cap \mathrm{lm}(\partial W_x^r)$ and $m \ne \mathrm{lm}(\partial_i pg)$ for any $i, g \in G$, and $p \in M_x$. For $m \in H$, let $a_m \in W_x^r$ be any reduced form of some $x^\gamma g - \mathrm{lc}(g)\partial^\beta x^{\alpha+\gamma}e_j$, where $g \in G$, $\mathrm{lm}(g) = x^\alpha \partial^\beta e_j$, and $m = \mathrm{lm}(x^\gamma g)$. Then

$$E_{\leq \eta} = \sum_{m \in H} \mathbb{K}a_m. \qquad (5)$$

*Proof.* Note that for each $m \in H$, the corresponding $\beta$ is nonzero. Indeed, by definition, $m \in \mathrm{lm}(\partial W_x^r)$ so there is some $\partial_i$ such that $m = \mathrm{lm}(\partial_i m')$ for another monomial $m'$. Moreover, $m = \mathrm{lm}(x^\gamma g)$, so $\mathrm{lm}(g)$ also contains $\partial_i$. In particular, the term $\mathrm{lc}(g)\partial^\beta x^{\alpha+\gamma}e_j$ has the form $\partial_i w$. Therefore, Theorem 6 applies and $E_{\leq m} = E_{\leq m} + \mathbb{K}a_m$. For a monomial $m$ not in $H$, either Theorem 6 or Lemma 7 shows that $E_{\leq m} = E_{\leq m}$. Then the statement follows by well-founded induction on $\eta$. $\square$

To turn Corollary 8 into an algorithm, we introduce a finiteness property of the monomial order $\leq$.

**Hypothesis 9.** For any two monomials $\gamma$ and $\eta$ of $M_{x,r}$, the set of $\alpha$ for which $x^\alpha \gamma \le \eta$ is finite.

This hypothesis is always satisfied by orders graded by total degree, because a monomial $\eta$ has a finite number of predecessors in $M_{x,r}$. It is also satisfied by orders eliminating $x$, in the sense that

$$\alpha' - \alpha \in \mathbb{N}^n \setminus \{0\} \Rightarrow x^\alpha \partial_x^\beta e_i \le x^{\alpha'} \partial_x^{\beta'} e_i, \qquad (6)$$

as long as the set of $\alpha$ for which $x^\alpha \le x_i$ is finite for each $i \in \{1, \dots, n\}$. For example, this contains “elimination orders” [25] or “block orders” [3] that first order by total degree in $x$, but not a lexicographical order that has $x_1 > x_2 > \partial_1 > \partial_2$.

**Theorem 10.** *Under Hypothesis 9 Algorithm 2 is correct and terminates.*

*Proof.* Termination is obvious since the set on line 2 is finite, by hypothesis. For the correction, we observe that $H$ computed in the algorithm is the same as the set $H$ described in Corollary 8. $\square$

**Definition 11.** Let $B_\eta$ be an echelon form of the generating family returned by Algorithm 2 on input $\eta$. We define a reduction $[\cdot]_\eta$ from $W_x^r$ into itself by

$$[a]_\eta = \mathrm{Reduce}([a], B_\eta)$$

where $[\cdot]$ is the map defined by Algorithm 1 and $\mathrm{Reduce}(., B_\eta)$ is the reduction algorithm by the echelon form $B_\eta$.
**Algorithm 2** Computation of $E_{\leq \eta}$

**Input:**

• a Gröbner basis $G$ of $S$

• $\eta \in M_{\mathbf{x},r}$

**Output:**

• a generating family of the $\mathbb{K}$-vector space $E_{\leq \eta}$

1. $G' \leftarrow \{g \in G \mid \mathrm{lm}(g) \text{ involves some } \partial_i\}$

2. $H \leftarrow \{\mathrm{lm}(\mathbf{x}^{\gamma}g) \mid \gamma \in \mathbb{N}^n, g \in G', \text{ and } \mathrm{lm}(\mathbf{x}^{\gamma}g) \leq \eta\} \# \text{ finite by Hypothesis 9}$

3. $H \leftarrow H \setminus \{\mathrm{lm}(\partial_i mg) \mid m \in M_{\mathbf{x}}, g \in G, 1 \leq i \leq n\}$

4. $B \leftarrow \emptyset$

5. **for** $m \in H$

6.   pick $g \in G'$ and $\gamma$ such that $m = \mathrm{lm}(\mathbf{x}^{\gamma}g)$

7.   $\mathbf{x}^{\alpha} \partial^{\beta} e_i \leftarrow \mathrm{lm}(g) \quad \# \text{ by construction } \beta \neq 0$

8.   $B \leftarrow B \cup \{[\mathbf{x}^{\gamma}g - \mathrm{lc}(g) \partial^{\beta} \mathbf{x}^{\alpha} + \gamma e_i]\}$

9. **return** $B$

**Proposition 12.** The map $[\cdot]_{\eta}$ is $\mathbb{K}$-linear.

*Proof.* This follows from Proposition 5 and the $\mathbb{K}$-linearity of Reduce($\cdot$, $B_{\eta}$). $\square$

**Theorem 13.** For any $a \in S + \partial W_{\mathbf{x}}^r$ there exists $\eta \in M_{\mathbf{x},r}$ such that for all $\eta' \geq \eta$, the remainder $[a]_{\eta'}$ is zero.

*Proof.* The element $[a]$ is congruent to $a$ modulo $S + \partial W_{\mathbf{x}}^r$, so it is in $S + \partial W_{\mathbf{x}}^r$, like $a$ itself. Moreover, it is irreducible, and so it is in $E$ by the definition of $E$. Because of the equality $E = \bigcup_{\eta \in M_{\mathbf{x},r}} E_{\leq \eta}$, there exists $\eta$ such that $[a] \in E_{\leq \eta}$ and thus $\mathrm{Reduce}([a], B_{\eta}) = 0$. For $\eta' \geq \eta$, the vector space $\mathrm{Span}_{\mathbb{K}}(B_{\eta})$ is included in $\mathrm{Span}_{\mathbb{K}}(B_{\eta'})$, so $[a]_{\eta'} = \mathrm{Reduce}([a], B_{\eta'}) = \mathrm{Reduce}(\mathrm{Reduce}([a], B_{\eta}), B_{\eta'}) = \mathrm{Reduce}(0, B_{\eta'}) = 0$. $\square$

**Definition 14.** The normal form of an element $a \in W_{\mathbf{x}}^r$ modulo $S + \partial W_{\mathbf{x}}^r$ is the unique element $a' \in W_{\mathbf{x}}^r$ such that $a \equiv a' \pmod{S + \partial W_{\mathbf{x}}^r}$ and no monomial of $a'$ is the leading monomial of an element of $S + \partial W_{\mathbf{x}}^r$.

**Corollary 15.** For any $a \in W_{\mathbf{x}}^r$, there exists $\eta \in M_{\mathbf{x},r}$ such that for all $\eta' \geq \eta$, the remainder $[a]_{\eta'}$ is the normal form of a modulo $S + \partial W_{\mathbf{x}}^r$.

*Proof.* Let $a'$ be the normal form of $a$. Let $\eta$ such that $[a - a']_{\eta} = 0$, given by Theorem 13. By definition, $[\cdot]_{\eta}$ replaces monomials by smaller ones, but only if this is possible, so that we have $[a']_{\eta} = a'$. By linearity of $[\cdot]_{\eta}$, we obtain $[a]_{\eta} = a'$. $\square$
<caption>**Algorithm 3** Computation of a confinement</caption>

**Input:**
* a module $S \subseteq W_x^r$ given by a Gröbner basis
* $a \in W_x^r$
* a $W_x$-linear map $L : W_x^r \to W_x^r$ given by a $r \times r$ matrix
* $\rho \in \mathbb{N}$

**Output:**
* an effective confinement for $a$ and $L$ modulo $S + \partial W_x^r$

1 &nbsp;&nbsp;&nbsp;&nbsp; $s \leftarrow \rho$
2 &nbsp;&nbsp;&nbsp;&nbsp; $\eta \leftarrow$ the largest monomial of degree $s$
3 &nbsp;&nbsp;&nbsp;&nbsp; $B \leftarrow \emptyset$
4 &nbsp;&nbsp;&nbsp;&nbsp; $Q \leftarrow \text{supp}([a]_\eta)$
5 &nbsp;&nbsp;&nbsp;&nbsp; **while** $Q \setminus B \neq \emptyset$
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $m \leftarrow$ an element of $Q \setminus B$
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $\text{deg } m > s - \rho$
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $s \leftarrow s + 1$
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **goto** line 2
10 &nbsp;&nbsp;&nbsp;&nbsp; $Q \leftarrow Q \cup \text{supp}([L(m)]_\eta)$
11 &nbsp;&nbsp;&nbsp;&nbsp; $B \leftarrow B \cup \{m\}$
12 &nbsp;&nbsp;&nbsp;&nbsp; **return** $(\eta, B)$

## 3.3 Confinement

“Computing” in the quotient $M/\partial M \simeq W_x^r/(S + \partial W_x^r)$ can take on several forms, with various levels of potency. In the strongest interpretation, we want to compute a basis of the quotient, as a $\mathbb{K}$-linear space, and we want to be able to compute normal forms in $W_x^r$ modulo $S + \partial W_x^r$. In a weaker sense, we merely want to be able to capture the finiteness of the quotient space, without ensuring the linear independence of a finite generating set or even producing it explicitly. In view of our needs for integration algorithms in the next sections, there is an even weaker sense: given $a \in W_x^r$ (which will designate a function to be integrated) and a $W_x$-linear map from $W_x^r$ to itself (which will be related to taking derivatives with respect to a parameter $t \in \mathbb{K}$), we need to testify the finite-dimensionality of the span over $\mathbb{K}$ of the orbit $\{L^i(a) \mid i \in \mathbb{N}\}$ modulo $S + \partial W_x^r$. In this section, we show that the reductions $[\cdot]_\eta$ can be used to find, for any $a$ and $L$, a finite set $B$ that witnesses this finite-dimensionality.

**Definition 16.** An *effective confinement* for $a \in W_x^r$ and a $W_x$-linear map $L : W_x^r \to W_x^r$ is a pair $(\eta, B)$ consisting of a monomial $\eta$ and of a finite subset $B \subseteq M_{x,r}$, and satisfying:

1. the support of $[a]_\eta$ is included in $B$;
2. the support of $[L(m)]_\eta$ is included in $B$ for any $m \in B$.

An effective confinement is free if the elements of $B$ are $\mathbb{K}$-linearly independent modulo $S + \partial W_x^r$.

**Theorem 17.** *Algorithm 3 is correct. It terminates if $M/\partial M$ is finite-dimensional, for example if $M$ is holonomic. Moreover, if the input parameter $\rho$ is large enough, then Algorithm 3 outputs a free effective confinement.*

*Proof.* We first address correctness. Consider the sets $B$ and $Q$ after any iteration of the while loop. By construction, we have $B \subseteq Q$, $\text{supp}([a]_\eta) \subseteq Q$, and $\text{supp}([L(m)]_\eta) \subseteq Q$ for any $m \in B$. If the halting condition $Q \setminus B \neq \emptyset$ of the while loop is reached, that is, equivalently, if $Q \subseteq B$ holds at the end of an iteration, then we have $B = Q$. In conclusion, the returned value $(\eta, B)$ is an effective confinement.

As for termination, let $C \subseteq M_x^r$ be the set of monomials that are normal forms modulo $S + \partial W_x^r$. As a consequence of Definition 14, the set of normal forms is the vector space $\text{Span}_{\mathbb{K}}(C)$. Moreover, a basis of the quotient space $M/\partial M$ is formed by the classes modulo $S + \partial W_x^r$ of all elements of $C$, so in particular, $C$ is finite by the hypothesis of finite dimension. By Corollary 15, there is therefore some $\eta_\infty$ such that for any $\eta \ge \eta_\infty$,

$$ \text{supp}([a]_\eta) \subseteq C \text{ and } \forall m \in C, \text{supp}([L(m)]_\eta) \subseteq C. \quad (7) $$

Since each iteration of the *while* loop treat a different monomial $m$, and since there are finitely many monomials of degree at most $s - \rho$, the *while* loop terminates. It terminates either because $Q \setminus B = \emptyset$, in which case the algorithm terminates, or because $Q$ contains an element of degree larger that $s - \rho$, in which case we increase $s$. So, either the algorithm terminates, or $s$ tends to $\infty$.

Assume $s \to \infty$. At some point, we will have $s \ge \text{deg} \eta_\infty$, so after line 2 is executed, we have the inequalities $\eta \ge \eta_\infty x_1^{s-\text{deg}(\eta_\infty)} \ge \eta_\infty$, because $\eta$ is the largest monomial of degree $s$ and by the definition of a monomial order. In this circumstances, the set $Q$ is a subset of $C$ at every iteration of the main loop, because of (7), and so is $B$ because of the invariant $B \subseteq Q$. Since $s \to \infty$, we also reach a point where $\rho + \text{deg} m \le s$ for all $m \in C$. After this point, $s$ is not increased anymore. This contradiction shows that the algorithm terminates.

If the input $\rho$ satisfies $\rho \ge \text{deg} \eta_\infty$, we have $s \ge \rho \ge \text{deg} \eta_\infty$, and by the same reasoning as in the previous paragraph, we have again $\eta \ge \eta_\infty$. This is so during the whole execution of the algorithm. Therefore, like in the preceding paragraph, we have $B \subseteq Q \subseteq C$ during the execution of the while loop. So, the output set $B$ is a subset of $C$, which is a free family modulo $S + \partial W_x^r$. $\square$

## 3.4 Comparison with the Griffiths–Dwork reduction

Let $f \in \mathbb{K}[\mathbf{x}]$ be a homogeneous polynomial and let $M$ be the $W_x$-module $\mathbb{K}[\mathbf{x}]e^f$, where $\partial_i$ acts by $\partial_i \cdot e^f = \frac{\partial f}{\partial x_i} e^f$. When $f$ defines a smooth variety, we can compute in $M/\partial M$ using the Griffiths–Dwork reduction [27, 28, 32]. This
is usually presented with rational functions in $K[x, f^{-1}]$ but the exponential formulation is equivalent (for example, see [26, 38]). The module $M$ admits the presentation

$$M \simeq \frac{W_x}{\sum_i W_x(\partial_i - f_i)},$$

where $f_i$ denotes the partial derivative $\frac{\partial f}{\partial x_i}$. (This presentation is what makes the exponential formulation easier in our setting. A presentation of the holonomic $W_x$-module $K[x, f^{-1}]$ is much harder to compute [44].)

We briefly present the Griffiths–Dwork reduction and observe that irreducible elements for the Griffiths–Dwork reduction are exactly the irreducible elements for our reduction $\to$.

**The Griffiths–Dwork reduction.** Let $\ll_0$ be a monomial order on $K[x]$, and, for this monomial ordering, let $G_0$ be the minimal Gröbner basis of the polynomial ideal $I = (f_1, \dots, f_n)$. Given a homogeneous polynomial $a \in K[x]$, we can compute the remainder $r$ of the multivariate division of $a$ by $G_0$ and the cofactors $b_1, \dots, b_n \in K[x]$ such that

$$a = r + \sum_{i=1}^{n} b_i f_i. \qquad (8)$$

By homogeneity, $\deg b_i = \deg a - \deg f + 1$ (unless $b_i = 0$). Then, the rule for the derivative of a product yields

$$ae^f = re^f - \underbrace{\sum_{i=1}^{n} \frac{\partial b_i}{\partial x_i} e^f}_{\text{degree } \deg a - \deg f} + \underbrace{\sum_{i=1}^{n} \partial_i \cdot (b_i e^f)}_{\in \partial M}. \qquad (9)$$

The last term is in $\partial M$, so we ignore it, and the second term has lower degree than $a$, so we can apply the same procedure recursively, which will terminate by induction on the degree. In the end, we obtain a reduced form $ae^f \equiv re^f \pmod{\partial M}$ where $r \in K[x]$ is irreducible with respect to $G_0$. These are the irreducible elements for the Griffiths–Dwork reduction.

This reductions is defined for any homogeneous polynomial $f$, but it enjoys special properties when $f_1, \dots, f_n$ do not have any non trivial common zero in an algebraic closure of $K$. Geometrically, this means that $f$ defines a smooth hypersurface in $\mathbb{P}^{n-1}(K)$. Griffiths [32] proved, under this smoothness assumption, that the reduced form of any $ae^f$ vanishes if and only if $ae^f \in \partial M$.

**Comparison with the reduction $\to$**. We consider the reduction rule $\to$ applied to the left ideal $S$ of $W_x$ generated by $\partial_i - f_i$, so that $K[x]e^f \simeq W_x/S$ (Section 3.1). Let $\ll$ be a monomial order on $W_x$ that eliminates $x$ (see (6)) and agrees with $\ll_0$ on $K[x]$. By following the steps of Buchberger’s algorithm, we observe that there is a Gröbner basis $G$ of $S$ in which each element is:
1. either an element of the form $r - \sum_{i=1}^{n} b_i \partial_i$, with $r \in G_0$, $b_i \in \mathbb{K}[\mathbf{x}]$ of degree $\deg(r) - \deg(f_i)$, and $r = \sum_i b_i f_i$,

2. or an element of $W_{\mathbf{x}}\partial_1 + \cdots + W_{\mathbf{x}}\partial_n$.

We will call such elements respectively of the first kind and of the second kind.

We now characterize irreducible elements in $W_{\mathbf{x}}$ with respect to $\to$. Let $a \in W_{\mathbf{x}}$ be an irreducible element. Since $a$ cannot be reduced with $\to_2$, it contains no $\partial_i$, so it is a polynomial. Since $a$ cannot be reduced with $\to_1$, no monomial in $a$ is divisible by the leading term of an element of $G$. By considering the elements of the first kind, we see that the monomials of $a$ are not divisible by the leading monomial of any element of $G_0$. So $a$ is irreducible with respect to $G_0$. The converse also holds: if $a \in \mathbb{K}[\mathbf{x}]$ is irreducible with respect to $G_0$, then $a$ is irreducible in $W_{\mathbf{x}}$ with respect to $S$. In this sense, we can regard Algorithm 1 as a generalization of the Griffiths-Dwork reduction.

As we observed, this reduction is not enough to compute in $M/\partial M$, since there may be nonzero irreducible elements in $S + \partial W_{\mathbf{x}}$. In the case of rational functions, Lairez [39] gave an algorithm to compute them efficiently. The algorithm that we have given in Section 3.2 behaves differently. In short, the algorithm in [39] would only consider elements of the second kind with degree 1 in the $\partial_i$, whereas we consider all elements of the second kind. On the one hand, this seems to give more reduction power, on the other hand the cost of computing them is higher. This indicates room for improvement in future work.

# 4 Creative Telescoping by Reduction

In the previous section, we obtained an algorithm for normalizing modulo derivatives in a holonomic $W_{\mathbf{x}}$-module. In this section, we introduce a parameter $t$ and differentiation with respect to $t$. It would be natural to work with a holonomic $W_{t,\mathbf{x}}$-module, but in view of the previous section, we need a finitely presented module over a Weyl algebra in the derivatives with respect to $\mathbf{x}$ only. This motivates the following context.

We consider the Weyl algebra $W_{\mathbf{x}}(t) = \mathbb{K}(t) \otimes_{\mathbb{K}} W_{\mathbf{x}}$ (which is just a Weyl algebra over the field $\mathbb{K}(t)$), and a holonomic $W_{\mathbf{x}}(t)$-module $M$ with a compatible derivation $\partial_t$, that is, a $\mathbb{K}$-linear map $\partial_t : M \to M$ such that for any $a \in W_{\mathbf{x}}(t)$ and any $m \in M$,

$$ \partial_t \cdot am = \frac{\partial a}{\partial t}m + a\partial_t \cdot m, $$

where $\frac{\partial a}{\partial t}$ is the coefficient-wise differentiation in $W_{\mathbf{x}}(t)$. In other words, $M$ is a $W_{t,\mathbf{x}}(t)$-module that is holonomic as a $W_{\mathbf{x}}(t)$-module.

It is also convenient to fix a finite presentation $W_{\mathbf{x}}(t)^r/S$ of $M$ and assume that there is a $W_{\mathbf{x}}(t)$-linear map $L : W_{\mathbf{x}}(t)^r \to W_{\mathbf{x}}(t)^r$ such that for any $a \in W_{\mathbf{x}}(t)^r$,

$$ \partial_t \cdot \mathrm{pr}_S(a) = \mathrm{pr}_S \left( \frac{\partial a}{\partial t} + L(a) \right), \quad (10) $$

where $\mathrm{pr}_S$ is the canonical map $W_{\mathbf{x}}(t)^r \to M$. In particular, note that $S$ is stable under $\frac{\partial}{\partial t} + L$. From the algorithmic point of view, we represent $M$ by
its finite presentation and the derivation $\partial_t$ by the the $r \times r$ matrix of the endomorphism $L$. We explain in Section 4.3 how to obtain this setting from a holonomic $W_{t,x}$-module.

Using the algorithm of the previous section, we aim at describing an algorithm that performs integration with respect to $x_1, \dots, x_n$, in the following sense. Given $f \in M$, we want to compute a nonzero operator $P(t, \partial_t) \in W_t(t)$ such that

$$P(t, \partial_t) \cdot f \in \partial M \qquad (11)$$

with the motivation that $P(t, \partial_t)$ is then an annihilating operator of the integral of $f$ with respect to $x_1, \dots, x_n$. The principle of integration by reduction is described in Section 4.1 and an algorithm is presented in Section 4.2.

## 4.1 Integration by reduction

We utilize the family of reductions $[\cdot]_\eta$ defined in Section 3.2. Let $f$ be an element of $W_x(t)^r$, let $\eta$ be some monomial in $M_{x,r}$ and let $(g_i)_{i \ge 0}$ be the sequence in $W_x(t)^r$ defined by

$$g_0 = [f]_\eta \quad \text{and} \quad g_{i+1} = \frac{\partial g_i}{\partial t} + [L(g_i)]_\eta \quad \text{for all } i \ge 0. \qquad (12)$$

As usual with integration-by-reduction algorithms, we relate the dependency relations between the reduced forms $g_i$ to the operators $P \in W_t(t)$ such that $P \cdot \text{pr}_S(f) \in \partial M$, which, as is traditional, we call telescopes for $f$.

**Lemma 18.** For any $i \ge 0$,

$$\text{pr}_S(g_i) \equiv \partial_t^i \cdot \text{pr}_S(f) \pmod{\partial M}. \qquad (13)$$

*Proof.* For $i=0$, this means that $[f]_\eta \equiv f \pmod{S + \partial W_x(t)^r}$, which holds by construction of $[\cdot]_\eta$. By property of $[\cdot]_\eta$, again, there is some $s_i \in S$ and some $\Delta_i = \sum_j \partial_j a_{i,j} \in \partial W_x(t)^r$ such that

$$[L(g_i)]_\eta = L(g_i) + s_i + \Delta_i.$$

Therefore, using the $W_x(t)$-linearity of $\text{pr}_S$ and $\text{pr}_S(s_i) = 0$, we obtain

$$\begin{align*}
\text{pr}_S(g_{i+1}) &= \text{pr}_S\left(\frac{\partial g_i}{\partial t} + L(g_i)\right) + \text{pr}_S(\Delta_i) \\
&= \partial_t \cdot \text{pr}_S(g_i) + \sum_j \partial_j \text{pr}_S(a_{i,j}), \quad \text{using (10)}, \\
&\equiv \partial_t \cdot \text{pr}_S(g_i) \pmod{\partial M}.
\end{align*}$$

The claim follows by induction on $i$, using that $\partial_t \cdot \partial M \subset \partial M$, since $\partial_t$ commutes with $\partial$.

**Lemma 19.** Let $P = \sum_{i=0}^N c_i \partial_t^i \in W_t(t)$.

1. If $c_0 g_0 + \dots + c_N g_N = 0$, then $P \cdot \text{pr}_S(f) \in \partial M$.
2. If $P \cdot \mathrm{pr}_S(f) \in \partial M$, then $c_0g_0 + \cdots + c_Ng_N \in S + \partial W_x(t)^r$.

3. If $\eta$ is large enough and if $P \cdot \mathrm{pr}_S(f) \in \partial M$, then $c_0g_0 + \cdots + c_Ng_N = 0$. (Note that "large enough" depends on $f$ and $P$.)

*Proof.* The first assertion follows directly from Lemma 18. Conversely, assume that $P \cdot \mathrm{pr}_S(f) \in \partial M$. This implies, again by Lemma 18, that

$$ \sum_{i=0}^{N} c_i g_i \in S + \partial W_x(t)^r, $$

proving the second assertion. Now, we observe that, if in addition $\eta$ is large enough, the $g_i$ are normal forms modulo $S + \partial W_x(t)^r$. Indeed, by Corollary 15, the $[L(g_i)]_\eta$ are normal forms; and since being a normal form is a condition on the monomial support, it is stable under coefficient-wise differentiation, so the $\frac{\partial g_i}{\partial t}$ are normal forms, by induction on $i$. So the linear combination $\sum_i c_i g_i$ is also a normal form modulo $S + \partial W_x(t)^r$. This implies $\sum_i c_i g_i = 0$. $\square$

## 4.2 An algorithm for integrating by reduction

To turn Lemma 19 into an algorithm to compute a telescoper, it only remains to find a suitable $\eta$. We use the idea of confinement (Section 3.3). Using Algorithm 3, we can compute an effective confinement for $f$ and $L$. Recall that this is a monomial $\eta$ and a finite set $B$ of monomials in $M_{x,r}$ such that $\mathrm{supp}([f]_\eta) \subseteq B$ and $\mathrm{supp}([L(b)]_\eta) \subseteq B$ for all $b \in B$. The following statement explains that reduced forms of successive derivatives with respect to $t$ therefore lie in the finite-dimensional vector-space $\mathrm{Span}_{\mathbb{K}(t)}(B)$.

**Lemma 20.** Let $(\eta, B)$ be an effective confinement for $f$ and $L$. Let $(g_i)_{i \ge 0}$ be the sequence defined by (12). Then, for all $i \ge 0$, $g_i \in \mathrm{Span}_{\mathbb{K}(t)}(B)$.

*Proof.* By definition of an effective confinement, $[f]_\eta \in \mathrm{Span}_{\mathbb{K}(t)}(B)$, and the space $\mathrm{Span}_{\mathbb{K}(t)}(B)$ is stable under $[L(\cdot)]_\eta$. Moreover, $\mathrm{Span}_{\mathbb{K}(t)}(B)$ is stable under $\frac{\partial}{\partial t}$. So the claim follows from the definition of $g_i$. $\square$

Algorithm 3 and Lemma 20 combine into Algorithm 4, whose main properties are provided in the following theorem.

**Theorem 21.** *Algorithm 4 is correct and terminates. Moreover, if $\rho$ is large enough, then it outputs a minimal telescoper for the input.*

*Proof.* Correctness follows from Lemma 19. As to termination, it follows from Lemma 20: because the set $B$ is finite, the infinite family of elements $g_0, g_1, \dots$ is linearly dependent, so the main loop terminates for some $N$ less than or equal to the cardinality of $B$.

As for the minimality, it is clear that the algorithm outputs a non-trivial relation $c_0g_0 + \cdots + c_Ng_N = 0$ with minimal possible $N$ among those available
**Algorithm 4** Integration using reductions

**Input:**

* a holonomic module $W_x(t)^r/S$

* a derivation map $\partial_t : W_x(t)^r/S \to W_x(t)^r/S$ given by the matrix of an endomorphism $L$, as in (10).

* an element $f \in W_x(t)^r$

* an integer $\rho \ge 0$

**Output:**

* $P = c_0 + \dots + c_N \partial_t^N$ such that $c_i \in \mathbb{K}(t)$, $c_N \ne 0$ and $P \cdot \mathrm{pr}_S(f) \in \partial M$

1. $(\eta, B) \leftarrow$ an effective confinement obtained from $(S, f, L, \rho)$ by Algorithm 3

2. $g_0 \leftarrow [f]_\eta$

3. $N \leftarrow 0$

4. **while** $g_0, \dots, g_N$ are linearly independent over $k(t)$

5. $g_{N+1} \leftarrow \frac{\partial g}{\partial t} + [L(g)]_\eta$

6. $N \leftarrow N + 1$

7. **return** $c_0 + \dots + c_N \partial_t^N$ s.t. $c_0 g_0 + \dots + c_N g_N = 0$, $c_i \in \mathbb{K}(t)$ and $c_N \ne 0$.

for the sequence $(g_i)_{i \ge 0}$. Besides, consider any telescoper $P = c_0 + \dots + c_\Omega \partial_t^\Omega$. By point 2 of Lemma 19, we have

$$c_0 g_0 + \dots + c_\Omega g_\Omega \in S + \partial W_x(t)^r. \quad (14)$$

Assume that $\rho$ is large enough, in the sense of Theorem 17, so that the confinement is free, meaning that the elements of $B$ are independent modulo $S + \partial W_x^r$. The linear combination in (14) is a linear combination of elements of $B$, so it must be zero: $c_0 g_0 + \dots + c_\Omega g_\Omega = 0$. So Algorithm 4 will output a relation for an $N$ that is at most the minimal order of telescopes. $\square$

**Remark 22.** Algorithm 4 can be modified to compute a system of linear differential equations satisfied by an integral depending on multiple parameters $t_1, \dots, t_p$. These parameters can also be associated to other Ore operators [19] than the differentiation provided they define a map on $W_x^r(t_1, \dots, t_p)$.

## 4.3 Scalar extension

Let $P$ be a holonomic $W_{t,x}$-module and let $M = \mathbb{K}(t) \otimes_{\mathbb{K}[t]} P$. This space $M$ is a $W_x(t)$-module in a natural way. Moreover, we can define a derivation $\partial_t$ by

$$\partial_t \cdot (a \otimes m) = \frac{\partial a}{\partial t} \otimes m + a \otimes (\partial_t \cdot m).$$

This derivation commutes with the action of the $\partial_i$.

In this section, we aim to compute $M$ from $P$, so that we can apply the integration algorithm (Section 4.2). Let us first make what we mean explicit.
We assume that $P$ is given by a finite presentation, that is, $P = W_{t,x}^s/J$ for some $s \ge 0$ and some submodule $J \subseteq W_{t,x}^s$ given by a finite set of generators. Computing $M$ means computing a finite presentation $M \simeq W_x(t)^r/S$, and a $W_x(t)$-linear map $L: W_x(t)^r \to W_x(t)^r$ such that (10) holds.

It is not obvious that such a finite presentation exists because $M$ does not have any obvious finite set of generators. However, this existence is implied by the holonomy of $M$. Here, we give a proof based on restriction of D-modules.

**Lemma 23.** If $P$ is a $W_{t,x}$-holonomic module, then $M = \mathbb{K}(t) \otimes_{\mathbb{K}[t]} P$ is a $W_x(t)$-holonomic module.

The statement is similar in nature to the well-known statement that “holonomic implies D-finite”.

*Proof.* Let $\xi$ be a new variable. Consider the field $\mathbb{L} = \mathbb{K}(\xi)$. Introduce the $(1+n)$th Weyl algebra with coefficients in $\mathbb{L}$, which we denote $W_{t,x}(\xi)$. Consider as well the $W_{t,x}(\xi)$-module $P' = \mathbb{L} \otimes_{\mathbb{K}} P$. This scalar extension of the base field this preserves holonomy. So $P'$ is holonomic. Consider now the embedding map $F: \mathbb{L}^n \to \mathbb{L}^{1+n}$ defined by

$$ (x_1, \dots, x_n) \mapsto (\xi, x_1, \dots, x_n), $$

and the inverse image $F^*P'$, that is, the restriction of $P'$ at $t = \xi$. We just need to know that $F^*P'$ is:

* a $W_x(\xi)$-module, which is by [24, construction of §14.1 and §14.2],

* holonomic as a $W_x(\xi)$-module, which is by [24, Theorem 18.1.4],

* isomorphic to $P'/(t-\xi)P'$, which is obtained as a suitable variant of [24, §15.1], by making $Y = t$ in that reference before specializing at $\xi$ instead of 0.

In particular, we have

$$ F^*P' \simeq P'/(t-\xi)P' \simeq \mathbb{L}[t]/(t-\xi) \otimes_{\mathbb{L}[t]} P'. $$

Next, we check that

$$ \begin{align*} P' &= \mathbb{L} \otimes_{\mathbb{K}} P, & \text{by definition} \\ &\simeq \mathbb{L} \otimes_{\mathbb{K}} (\mathbb{K}[t] \otimes_{\mathbb{K}[t]} P), & \text{because } P \text{ is a } \mathbb{K}[t] \text{ module} \\ &\simeq \mathbb{L}[t] \otimes_{\mathbb{K}[t]} P, & \text{by associativity of } \otimes, \end{align*} $$

and therefore, using associativity of $\otimes$ again,

$$ F^*P' \simeq \mathbb{L}[t]/(t-\xi) \otimes_{\mathbb{K}[t]} P. $$

Finally, we observe the isomorphism $\mathbb{K}(t) \simeq \mathbb{L}[t]/(t-\xi)$ as $\mathbb{K}[t]$-algebras under the map $f(t) \mapsto f(\xi)$, so we obtain $F^*P' \simeq \mathbb{K}(t) \otimes_{\mathbb{K}[t]} P = M$, by definition of $M$. Since $F^*P'$ is holonomic, this gives the claim. $\square$
We now describe an algorithm for computing $M$. Recall that $P = W_{t,\mathbf{x}}^s/J$, so $M \simeq W_{t,\mathbf{x}}(t)^s/J(t)$ where $J(t)$ is the submodule $K(t) \otimes_{K[t]} J$ of $W_{t,\mathbf{x}}(t)^s$ generated by $J$. We can compute normal forms in $M$ using Gröbner bases in $W_{t,\mathbf{x}}(t)^s$ after fixing a monomial order on all monomials $\mathbf{x}^\alpha \partial^\beta \partial_t^k e_i$ [e.g., 19]. We choose a monomial order that eliminates $\partial_t$, that is, any monomial order such that

$$k < k' \Rightarrow \mathbf{x}^{\alpha} \partial_{\mathbf{x}}^{\beta} \partial_{t}^{k} e_{i} < \mathbf{x}^{\alpha'} \partial_{\mathbf{x}}^{\beta'} \partial_{t}^{k'} e_{i'}.$$

Let $G$ denote a Gröbner basis of $J(t)$ for such an elimination order.

As a $W_{\mathbf{x}}(t)$-module, $W_{t,\mathbf{x}}(t)^s$ is generated by the set

$$\{\partial_t^i e_j \mid i \ge 0, 1 \le j \le s\}.$$

So $M$ is generated by the image of this set. This is an infinite family, but, since $M$ is $W_{\mathbf{x}}(t)$-holonomic, $M$ is actually a Noetherian $W_{\mathbf{x}}(t)$-module, finitely generated in particular. To describe $M$, we need to find an explicit finite generating set and the module of relations between the generators.

For $a \in W_{t,\mathbf{x}}(t)^s$, let $\mathrm{ind}(a)$ denote the degree of $a$ with respect to $\partial_t$, which we will call the index of $a$. In other words, this is the smallest integer $k \ge 0$ such that $a$ is in the sub-$W_{\mathbf{x}}(t)$-module generated by

$$B_k = \{\partial_t^i e_j \mid 0 \le i \le k, 1 \le j \le s\}.$$

Moreover, for $a \in W_{t,\mathbf{x}}(t)^s$, let $\mathrm{ind}_{J(t)}(a)$ denote

$$\mathrm{ind}_{J(t)}(a) = \min \{\mathrm{ind}(b) \mid b \in W_{t,\mathbf{x}}(t)^s \text{ and } a \equiv b \pmod{J(t)}\}. \quad (15)$$

Given $a$, we can compute $\mathrm{ind}_{J(t)}(a)$ using the Gröbner basis $G$:

$$\mathrm{ind}_{J(t)}(a) = \mathrm{ind}(\mathrm{LRem}(a, G)).$$

Indeed: we have $a \equiv \mathrm{LRem}(b, G)$ if $a \equiv b \pmod{J(t)}$ and the elimination property shows $\mathrm{ind}(\mathrm{LRem}(b, G)) \le \mathrm{ind}(b)$, so that $\mathrm{ind}(b)$ can be replaced with $\mathrm{ind}(\mathrm{LRem}(b, G))$ in (15); then the Gröbner basis property shows $\mathrm{LRem}(a, G) = \mathrm{LRem}(b, G)$ if $a \equiv b \pmod{J(t)}$.

**Lemma 24.** There is $\ell \ge 0$ such that $\mathrm{ind}_{J(t)}(\partial_t^{\ell+1} e_i) \le \ell$ for any $1 \le i \le s$. Moreover, for any such $\ell$:

1. $M$ is generated as a $W_{\mathbf{x}}(t)$-module by the image in it of $B_\ell$,

2. $\mathrm{ind}_{J(t)}(a) \le \ell$ for any $a \in W_{t,\mathbf{x}}(t)^s$.

*Proof.* Since $M$ is Noetherian, the increasing sequence of the $W_{\mathbf{x}}(t)$-modules $S_k$ generated by the images of the $B_k$ in $M$ is stationary: there exists $\ell \ge 0$ for which the $W_{\mathbf{x}}(t)$-module $S_\ell$ contains all the $S_k$ for $k \ge 0$, and is therefore equal to $M$. For such an integer $\ell$, any $k > \ell$, and any $j$, the image of $\partial_t^k e_j$ in $M$ is in $S_k$, therefore in $S_\ell = M$. Consequently, there exist coefficients $c_{h,i} \in W_{\mathbf{x}}(t)$ satisfying

$$\partial_t^k e_j \equiv \sum_{h \le \ell, i} c_{h,i} \partial_t^h e_i \pmod{J(t)}.$$
By the definition (15), ind$_{J(t)}$(∂_t^k e_j) is less than or equal to the index of the right-hand side, which by construction is less than or equal to ℓ. We obtain that ℓ is a uniform bound on all ind$_{J(t)}$(∂_t^k e_j). This proves in particular the first part of the statement, on the existence of ℓ. For the second part, we fix such an ℓ. We have already proved $M = S_ℓ$, which is the first itemized statement. We have also already proved, for any $k > ℓ$ and any $j$, the existence of some $r_{k,j}$ of index at most $ℓ$ such that $∂_t^k e_j ≡ r_{k,j}$ (mod $J(t)$). This also holds by the definition of the index for $k ≤ ℓ$. Now, any $a ∈ W_{t,x}^s$ writes in the form $\sum_{k,j} c_{k,j} ∂_t^k e_j$ for coefficients $c_{k,j} ∈ W_x(t)$. Taking linear combinations of congruences modulo the $W_x(t)$-module $J(t)$, we obtain $a ≡ \sum_{k,j} c_{k,j} r_{k,j}$ (mod $J(t)$), and as a consequence,

$$ \mathrm{ind}_{J(t)}(a) \le \mathrm{ind}\left(\sum_{k,j} c_{k,j} r_{k,j}\right) \le l, $$

where the first inequality is by (15) and the second by the definition of the index as a degree. We have proved the second itemized statement. □

An algorithm for computing the smallest $ℓ$ as in the statement above follows directly from (15), simply by testing increasing values of $ℓ$.

Now that we have a finite generating set for $M$, it remains to characterize the relations between the generators. To this end, for the rest of the section we fix $ℓ$ as provided by Lemma 24 and we let $J_ℓ$ be the sub-$W_x(t)$-module of $W_{t,x}(t)^s$ generated by

$$ \{\partial_t^k g \mid g \in G \text{ and } k + \mathrm{ind}(g) \le l\}. $$

It is, by construction, a submodule of $W_x(t)B_ℓ$.

**Lemma 25.** The inclusion $W_x(t)B_ℓ \to W_{t,x}(t)^s$ induces an isomorphism

$$ M \simeq \frac{W_x(t)B_\ell}{J_\ell}, $$

with inverse induced by the map $W_{t,x}(t)^s \to W_x(t)B_ℓ$ given by $a \mapsto \mathrm{LRem}(a, G)$.

*Proof.* First, $J_ℓ \subseteq J(t)$, so the inclusion $W_x(t)B_ℓ \to W_{t,x}(t)^s$ induces a morphism of $W_x(t)$-modules

$$ \phi : W_x(t)B_\ell/J_\ell \to W_{t,x}(t)^s/J(t). $$

Next, the $K(t)$-linear map $a \in W_{t,x}(t)^s \mapsto \mathrm{LRem}(a, G)$ has values in $W_x(t)B_ℓ$, because ind($\mathrm{LRem}(a, G)$) = ind$_{J(t)}(a) \le ℓ$, by Lemma 24. This map vanishes on $J(t)$, because $G$ is a Gröbner basis of $J(t)$, so it induces a $K(t)$-linear map

$$ \psi : W_{t,x}(t)^s/J(t) \to W_x(t)B_\ell/J_\ell. $$

The maps $\phi \circ \psi$ and $\psi \circ \phi$ are both induced by $a \mapsto \mathrm{LRem}(a, G)$. The first is the identity on $W_{t,x}(t)^s/J(t)$ because for all $a \in W_{t,x}^s$, $a \equiv \mathrm{LRem}(a, G) (\mathrm{mod} \ J(t))$ as a property of the Gröbner basis $G$. The second is the identity on $W_x(t)B_\ell/J_\ell$
because of the elimination property: the computation of LRem($a, G$) only involves multiples of $G$ of index at most ind($a$), which are all in $J_ℓ$, so that for all $a$ of index at most $ℓ$, $a \equiv \text{LRem}(a, G) (\text{mod } J_ℓ)$. This shows that $\phi$ is an isomorphism. $\square$

At this point we are able to define the wanted dimension $r$ and module $S$. In view of Lemma 25, set $r$ to $(ℓ + 1)s$, so as to have a trivial isomorphism $W_\mathbf{x}(t)B_ℓ \simeq W_\mathbf{x}(t)^r$. Call $S$ the image of the submodule $J_ℓ$ of $W_\mathbf{x}(t)B_ℓ$ under this isomorphism, so that, summarizing,

$$ \frac{W_{t,\mathbf{x}}^s}{J(t)} \simeq M \simeq \frac{W_{\mathbf{x}}(t)B_\ell}{J_\ell} \simeq \frac{W_{\mathbf{x}}(t)^r}{S}. $$

It remains to describe an endomorphism $L$ of $W_\mathbf{x}(t)^r$ such that

$$ \partial_t \cdot \text{pr}_S(a) = \text{pr}_S \left( \frac{\partial a}{\partial t} + L(a) \right), $$

for any $a \in W_\mathbf{x}(t)^r$. Introduce the canonical maps $\text{pr}_{J(t)}$ and $\text{pr}_S$ to the relevant quotients. Recall that $\partial_t$ is defined for any $h$ in $W_{t,\mathbf{x}}(t)^s/J(t)$ by left-multiplication by $\partial_t$:

$$ \partial_t \cdot \text{pr}_{J(t)}(h) = \text{pr}_{J(t)}(\partial_t h). $$

Therefore, the isomorphism of Lemma 25 transfers $\partial_t$ on $W_\mathbf{x}(t)B_\ell/J_\ell$ as

$$ \partial_t \cdot \text{pr}_{J_\ell}(a) = \text{pr}_{J_\ell}(\text{LRem}(\partial_t a, G)). \qquad (16) $$

Lastly, take $a \in W_\mathbf{x}(t)B_\ell$ and write it $a = \sum_{i=1}^n a_i e_i$. The Leibniz rule in $W_{t,\mathbf{x}}^s$ gives $\partial_t a = \frac{\partial a}{\partial t} + \sum_{i=1}^n a_i \partial_t e_i$. By linearity of LRem and since $W_\mathbf{x}(t)B_\ell$ is stable under the coefficient-wise differentiation $\frac{\partial}{\partial t}$, we obtain

$$ \begin{aligned} \text{LRem}(\partial_t a, G) &= \text{LRem}(\frac{\partial a}{\partial t}, G) + \text{LRem}\left(\sum_{i=1}^n a_i \partial_t e_i, G\right) \\ &= \frac{\partial a}{\partial t} + \text{LRem}\left(\sum_{i=1}^n a_i \partial_t e_i, G\right) + h \end{aligned} $$

for some $h \in J_\ell$, then, upon applying $\text{pr}_{J_\ell}$ and combining with (16),

$$ \partial_t \cdot \text{pr}_{J_\ell}(a) = \text{pr}_{J_\ell} \left( \frac{\partial a}{\partial t} + \text{LRem} \left( \sum_{i=1}^{n} a_i \partial_t e_i, G \right) \right). $$

So, the endomorphism $L$ we want is obtained by transferring the endomorphism

$$ \sum_{i=1}^{n} a_i e_i \mapsto \text{LRem} \left( \sum_{i=1}^{n} a_i \partial_t e_i, G \right) $$

of $W_\mathbf{x}(t)B_\ell$ to an endomorphism of $W_\mathbf{x}(t)^r$.
# 5 Implementation

This section outlines specific algorithmic and implementation choices made in our Julia implementation of the algorithms presented in this paper. We begin in Section 5.1 with a general presentation of our package. In Section 5.2 we present two ideas based on memoization and on the use of a tracer to efficiently compute the image of the reduction map $[\cdot]_\eta$ on a finite-dimensional space. Lastly, we present in Section 5.3 a modified version of Algorithm 4 that utilizes the evaluation/interpolation paradigm to avoid the growth of intermediate coefficients. This is of particular importance when reducing operators with coefficients in $K(t)$ using the reduction $[\cdot]_\eta$.

## 5.1 General comments

We have implemented our algorithms in Julia. This is available as the package `MultivariateCreativeTelescoping.jl`<fnref n="1" />. The package includes an implementation of Weyl algebras in which operators have a sparse representation by a pair of vectors, one for exponents and one for the corresponding coefficients. The currently supported coefficient fields are the field $\mathbb{Q}$ of rational numbers, the finite fields $\mathbb{F}_p$ with $p \le 2^{31}$, and extensions of those fields with symbolic parameters. When such parameters are present, the implementation interfaces with FLINT for defining and manipulating commutative polynomials, by means of the Julia packages `AbstractAlgebra.jl` and `Nemo.jl`. Our package also provides an implementation of non-commutative generalizations of algorithms for computing Gröbner bases: Buchberger's algorithm (see e.g. [2]), the F4 algorithm [29], and the F5 algorithm [47, 40]. The F4 implementation seems to be the most efficient one for grevlex orders while the F5 implementation seems to be the most efficient one for block and lexicographical orders.

## 5.2 Efficient computation of the reduction map $[\cdot]_\eta$ on a finite-dimensional vector space

We first explain how to compute efficiently the sequence $(g_i)_{i \ge 0}$ that is contained in a finite-dimensional vector space obtained by an effective confinement. We next discuss the use of a tracer to compute the vector space $E_{\le \eta}$. We finally provide comments on the dimension of $E_{\le \eta}$ for various examples.

**Computation of the sequence $(g_i)_i$ using memoization.** The confinement $(\eta, B)$ required by Algorithm 4 is constructed so that the sequence $(g_i)_i$ defined in (12) is contained in $\text{Span}_{K(t)}(B)$. Properties of the reduction $[\cdot]_\eta$ imply a refined formula: after decomposing $g_i$ as a sum $\sum_{m \in B} a_m m$ with $a_m \in K(t)$, $g_{i+1}$ can be obtained by

$$g_{i+1} = \sum_{m \in B} \left( \frac{\partial a_i}{\partial t} m + a_m [L(m)]_\eta \right). \qquad (17)$$

<footnote n="1">¹See https://hbrochet.github.io/MultivariateCreativeTelescoping.jl/.</footnote>
Early in the execution of Algorithm 4, when it computes the confinement $(\eta, B)$ by Algorithm 3, the image $[L(m)]_\eta$ of every monomial $m \in B$ has to be computed. We therefore choose to store these images in memory to allow for a more efficient computation of the sequence $(g_i)_i$ by using (17) at a later stage.

**The vector space $E_{\prec\eta}$.** In Algorithm 2, a generating set of the vector space $E_{\prec\eta}$ is computed by reducing for each $\eta' \prec \eta$ a term of the form $mg - \partial_i w$ satisfying $\eta' = \text{lm}(mg) = \text{lm}(\partial_i w)$. However, not all such terms contribute to an increase in the dimension of the space, as their reductions may be linearly dependent. Since such reductions are repeated multiple times for different primes $p$ and evaluation points of $t$ (see the next subsection), we use a tracer [51] during the computation for the first pair $(p, t)$ to record all $\eta'$ corresponding to a non-contributing term and skip the corresponding terms in subsequent computations. Assuming that the first pair $(p, t)$ is not unlucky, we know that all the skipped pairs would also lead to unnecessary elements if used in later computations.

**Dimension of the vector space $E_{\prec\eta}$.** By the finite dimensionality of $M/\partial M$, only finitely many monomials of $M_{x,r}$ are irreducible modulo $S + \partial W_x^r$. As a consequence, every monomial $m \in M_{x,r}$ except for the finitely many irreducible ones is either reducible by $G$ or by an echelon form of $E_{\prec\eta}$ for some $\eta$. If an infinite number of monomials is not reducible by $G$, as in Example 3, the dimension of $E_{\prec\eta}$ will tend to infinity when $\eta$ increases indefinitely, making the computation of $E_{\prec\eta}$ increasingly expensive. As a consequence, the computational cost of Algorithm 2 depends on the structure of the staircase formed by the leading monomials of $G$. We present two extreme scenarios: in one, $E_{\prec\eta}$ is equal to $\{0\}$ for any $\eta$ and in the other, no monomial of $K[x]$ is reducible by $G$. Naturally, intermediate cases also exist.

*Example 26.* Let $S$ be the left ideal of $W_x(t)$ generated by the Gröbner basis

$$ (t - 1)\underline{x}_1 - t\partial_1, \quad \underline{x}_2 - t. $$

Up to renaming variables, this is the left ideal used for the computation of the generating series of 2-regular graphs in Section 6. Every operator in this Gröbner basis has its leading monomial in $K[x]$, therefore Theorem 6 implies that $E_{\prec\eta}$ is $\{0\}$ for any $\eta$. In this very special case we obtain that the reduction [.] computes normal forms. That is, $[a] = 0$ if and only if $a \in S + \partial W_x^r$. We observed the same phenomenon with the ideals $S$ defined in Theorem 29 for $k$-regular graphs up to $k = 8$.

*Example 27.* Let $f(x, t) = 1 - (1 - x_1x_2)x_3 - tx_1x_2x_3(1 - x_1)(1 - x_2)(1 - x_3)$ and set the context for $n = 3$. The integral of $1/f$ is related to the generating function of Apéry numbers [5]. We were able to compute, by a method that we do not describe here, a Gröbner basis for the grevlex order of a $W_x(t)$-ideal $S$ included in $\text{Ann}(1/f)$ such that $W_x(t)/S$ is holonomic. This Gröbner basis contains 26 operators but all of their leading monomials contain a $\partial_i$. Hence,
no monomial of $K[x]$ is reducible the Gröbner basis of $S$. We observed the same phenomenon for every rational function that we tried.

Lastly, we present the simplest example on which our reduction $[\cdot]_\eta$ is inefficient.

*Example 28.* We continue Example 3, in which we set $S$ to $W_{x_1} \partial_1$ and $\prec \leq$ to the lexicographic order $\partial_1 \prec x_1$. We remarked that $E \subset K[x_1]$, and reciprocally the equality $(i+1)x_1^i = \partial_1 x_1^{i+1} - x_1^{i+1} \partial_1$ proves that $K[x_1] \subset E$. The reduction $[\cdot]$ does not see that elements of $K[x_1]$ are reducible by $S + \partial W_{x_1}$ and Algorithm 2 ends up calculating a basis of $K[x_1]_{\leq \eta}$.

## 5.3 Modular methods and evaluation/interpolation

In this subsection, we fix $K = \mathbb{Q}$. We first recall the principle of the evaluation/interpolation paradigm and then we present a modified version of Algorithm 4 that incorporates this paradigm.

**The principle.** Let $r_1, \dots, r_\ell \in \mathbb{Q}(t)$ be rational functions and let $F : \mathbb{Q}(t)^\ell \to \mathbb{Q}(t)$ be a function computable using only additions, multiplications, and inversions. The rational function $F(r_1, \dots, r_\ell)$ can be reconstructed as an element of $\mathbb{Q}(t)$ from its evaluations $F(r_1(a_i), \dots, r_\ell(a_i)) \bmod p_j$ at several points $a_i \in \mathbb{F}_{p_j}$ and for several prime integers $p_j$ in three steps:

1. for each $j$, reconstruct $F(r_1, \dots, r_\ell) \bmod p_j$ in $\mathbb{F}_{p_j}(t)$ by Cauchy interpolation [31, Chapter 5.8],

2. reconstruct $F(r_1, \dots, r_\ell) \bmod N$ in $\mathbb{F}_N(t)$ where $N = \prod_j p_j$ by the Chinese remainder theorem [31, Chapter 5.4],

3. lift the integer coefficients of $F(r_1, \dots, r_\ell) \bmod N$ in $\mathbb{Q}$ by rational reconstruction [31, Chapter 5.10].

The first (resp. third) step requires a bound on the degree in $t$ of the result (resp. on the size of its coefficients) to determine the number of evaluations needed to ensure correctness. Since we do not have access to such bounds, we rely instead on a probabilistic approach: after successfully obtaining a result with a certain number of evaluations in step 1 (resp. 3), we compute one additional evaluation and check consistency with the previously obtained result.

Finally, this method can fail if some bad evaluation points or some bad primes are chosen. However, these situations are very rare, provided the prime numbers are sufficiently large and the evaluation points are selected uniformly at random over $\mathbb{F}_p$.
<caption>**Algorithm 5** Integration with evaluation/interpolation</caption>

**Input:**
* a holonomic module $W_x(t)^r/S$
* a derivation map $\partial_t : W_x(t)^r/S \to W_x(t)^r/S$ given by the matrix of an endomorphism $L$, as in (10)
* an element $f \in W_x(t)^r$
* an integer $\rho \ge 0$

**Output:**
* $P = c_0 + \dots + c_N \partial_t^N$ such that $c_i \in \mathbb{Q}(t)$, $c_N \ne 0$ and $P \cdot \text{pr}_S(f) \in \partial M$

1 **for** large random prime numbers $p = p_1, p_2, \dots$
2 &nbsp;&nbsp;&nbsp; $\bar{L} \leftarrow$ the endomorphism of $\mathbb{F}_p \otimes_{\mathbb{Z}} W_x^r(t)$ induced by $L$ via the reduction modulo $p$
3 &nbsp;&nbsp;&nbsp; **for** random numbers $a = a_1, a_2, \dots$ in $\mathbb{F}_p$
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\tilde{S} \leftarrow$ image of $S$ under evaluation at $t = a$ and reduction modulo $p$
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\tilde{L} \leftarrow$ the endomorphism of $\mathbb{F}_p \otimes_{\mathbb{Z}} W_x^r$ induced by $\bar{L}$ via the same evaluation
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$(\eta, B) \leftarrow$ an effective confinement obtained from $(\tilde{S}, \tilde{f}, \tilde{L}, \rho)$ by Algorithm 3 over $\mathbb{K} = \mathbb{F}_p$
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\tilde{g}_0 \leftarrow [\tilde{f}]_\eta$ where the reduction is over $\mathbb{K} = \mathbb{F}_p$
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;store $\tilde{g}_0$ as well as the images $[\tilde{L}(m)]_\eta$ that have been computed at line 6 for all $m \in B$
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;interpolate the coefficients of $\tilde{g}_0$ and of each $[\tilde{L}(m)]_\eta$ by elements of $\mathbb{F}_p(t)$
10 &nbsp;&nbsp;&nbsp;$N \leftarrow 0$
11 &nbsp;&nbsp;&nbsp;**while** $\bar{g}_0, \dots, \bar{g}_N$ are linearly independent over $\mathbb{F}_p(t)$
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\bar{g}_{N+1} \leftarrow \frac{\partial \bar{g}_N}{\partial t} + [\bar{L}(\bar{g}_N)]_\eta$ where the reduction is over $\mathbb{K} = \mathbb{F}_p(t)$
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$N \leftarrow N + 1$
14 &nbsp;&nbsp;&nbsp;store coefficients of the minimal non-trivial relation $\bar{c}_0\bar{g}_0 + \dots + \bar{c}_N\bar{g}_N = 0$ for $\bar{c}_i \in \mathbb{F}_p(t)$
15 reconstruct the coefficients $c_0, \dots, c_N$ in $\mathbb{Q}(t)$ from their values in the $\mathbb{F}_{p_j}(t)$
16 **return** $c_0 + \dots + c_N \partial_t^N$
**Implementation.** The evaluation/interpolation scheme described above can- not be directly applied to Algorithm 4 as it involves not only additions, multipli- cations and inversions, but also differentiations. Indeed, recall that the sequence ( g$_{i}$ )$_{i}$ is defined for η ∈ M$_{x,r}$ by g$_{0}$ = [ f ]$_{η}$ and

$$
g_{i+1} = \frac{\partial g_i}{\partial t} + [L(g_i)]_\eta
$$

where ∂g/∂t denotes the coefficient-wise differentiation of g in the basis M$_{x,r}$ . This differentiation does not commute with the evaluation of t , which prevents the sequence ( g$_{i}$ )$_{i}$ from being computed by evaluation/interpolation.

Algorithm 5 computes the matrix of the K ( t )-linear map [ L ( . )]$_{η}$ using eval- uation/interpolation and uses it to find a linear relation among the elements of the sequence ( g$_{i}$ )$_{i}$. The reconstruction of the coefficients from F$_{p}$ ( t ) to Q ( t ) is delayed to the end of the algorithm. We adopt the following conventions: the projection of an element g ∈ W$_{x}$ ( t ) r in F$_{p}$ ⊗$_{Z}$ W r x by evaluation at a point a and reduction modulo p is stored in a variable ˜ g and its projection in F$_{p}$ ⊗$_{Z}$ W r x ( t ) by reduction modulo p is stored in a variable ¯ g .

6 Application to the computation of ODEs sat- isfied by k -regular graphs

In this section, we illustrate our new algorithm with computations on a fam- ily of multivariate integrals of combinatorial origin. We compute linear ODEs satisfied by various models of k -regular graphs and generalizations. A distin- guishing feature of these integration problems is that they cannot be solved by classical creative telescoping algorithms, which perform computations over the field of rational functions in all variables: because the objects to be integrated have polynomial torsion, they are not functions, and such calculations would erroneously result in a zero integral.

ODEs for k ≤ 5 were obtained 20 years ago by naive linear algebra and elim- ination by Euclidean divisions [23]. This has recently been extended to k ≤ 7 by a multivariate analog of the reduction-based algorithm [11] in which the re- duction is modulo the polynomial image of several differential operators. Here we use our new algorithm to achieve k = 8.

6.1 Statement of the integration problem

We first briefly introduce the problem and its solution by an integral represen- tation. We refer the reader to [23] for further motivation, history, and details. Given a fixed integer k ≥ 2, a k -regular graph is a graph whose vertices all have degree k , that is, all have exactly k neighbors. We are interested in the enumerative generating function

$$
R_k(t) = \sum_{n \ge 0} r_{k,n} \frac{t^n}{n!},
$$
where $r_{k,n}$ is the number of $k$-regular labeled graphs on $n$ vertices. It is well known that the generating function $R_k(t)$ satisfies a linear ODE with polynomial coefficients in $t$.

In this section, we show how such a differential equation can be obtained from Algorithm 4. To this end, we use the classical formulation of $R_k(t)$ as a scalar product of two exponentials,

$$R_k(t) = \langle e^f, e^{tg} \rangle,$$

where $f$ and $g$ are explicit polynomials in indeterminate $p_1, \dots, p_k$ [22, Algorithm 1, Step a.], and where the scalar product is a classic tool in the combinatorics of symmetric functions. The scalar product is first defined on monomials by

$$\langle p_1^{r_1} \cdots p_k^{r_k}, p_1^{s_1} \cdots p_k^{s_k} \rangle = z_r \delta_{\mathbf{r},\mathbf{s}} \quad \text{for} \quad z_r = r_1! 1^{r_1} r_2! 2^{r_2} \cdots r_k! k^{r_k}. \quad (18)$$

With $z_r$ indexed by exponents in $\mathbb{N}^k$, we depart from the equivalent indexing $z_\lambda$ by partitions $\lambda_1 \ge \lambda_2 \ge \dots$ that is used classically as well as in [22]. The scalar product is then extended by bilinearity to left arguments in $\mathbb{Q}[[\mathbf{p}]]$ and right arguments in $\mathbb{Q}[\mathbf{p}]((t))$, making the scalar product live in $\mathbb{Q}((t))$. Note that, because in the symmetric-function theory the power function $p_i$ denotes the sum $x_1^i + x_2^i + \dots$, we use the more traditional $\mathbf{p}$ instead of $\mathbf{x}$ for the variables, in accordance with existing literature. Also, by $\mathbb{Q}[\mathbf{p}]((t))$ we mean the ring of formal sums with coefficients in $\mathbb{Q}[\mathbf{p}]$, with finitely many exponents towards $-\infty$ and potentially infinitely many towards $+\infty$. This is not a field. The use of Algorithm 4 is justified by the following statement.

**Theorem 29.** Let $f, g \in \mathbb{Q}[\mathbf{p}]$. Let $S$ be the left ideal of $W_{\mathbf{p}}(t)$ generated by

$$p_i - t \frac{\partial g}{\partial p_i}(u_1, \dots, u_k),$$

for $1 \le i \le k$, where $u_i = i(\frac{\partial f}{\partial p_i} - \partial_i)$. Then, $W_{\mathbf{p}}(t)/S$ is holonomic as a $W_{\mathbf{p}}(t)$-module. Write $\mathrm{pr}_S$ for the canonical projection $\mathrm{pr}_S: W_{\mathbf{p}}(t) \to W_{\mathbf{p}}(t)/S$. Then, $W_{\mathbf{p}}(t)/S$ can be endowed with a derivation $\partial_t$ commuting with $\mathbf{p}$ and $\partial_{\mathbf{p}}$ satisfying

$$\partial_t \cdot \mathrm{pr}_S(a) = \mathrm{pr}_S \left( \frac{\partial a}{\partial t} + ag(u_1, \dots, u_k) \right). \quad (19)$$

On input the module $W_{\mathbf{p}}(t)/S$ (for $r=1$), the derivation (19) (for the implied endomorphism $L: a \mapsto ag(u)$), the element $f = 1 \in W_{\mathbf{p}}(t)$, and any $\rho \ge 0$, Algorithm 4 outputs a nonzero differential operator $P(t, \partial_t)$ such that

$$P(t, \partial_t) \cdot \langle e^f, e^{tg} \rangle = 0.$$

The rest of the section is a proof of this statement. In particular, we fix the ideal $S$ as in theorem, and the holonomy of $W_{\mathbf{p}}(t)/S$ will be proven as Lemma 33, the existence of the derivation $\partial_t$ will be proven as Lemma 31, and the correctness Theorem 21 will prove the final result on the operator $P$.
We begin with a few preliminary definitions. First, given $h = h(p_1, \dots, p_k) \in \mathbb{Q}[\mathbf{p}]$, we write $\tilde{h}$ for $h(1p_1, 2p_2, \dots, kp_k)$. Second, we define a formal Laplace transform on monomials by

$$ \mathcal{L}(p_1^{r_1} \dots p_k^{r_k}) = r_1! p_1^{-r_1-1} \dots r_k! p_k^{-r_k-1} \quad (20) $$

and extend it by linearity into a map from $\mathbb{Q}[\mathbf{p}]$ to $\mathbb{Q}[\mathbf{p}^{-1}]$. Third, we introduce the ring

$$ \mathbb{Q}\langle\langle\mathbf{p}\rangle\rangle := \mathbb{Q}[[\mathbf{p}]][\mathbf{p}^{-1}] = \bigcup_{\ell \ge 0} (p_1 \cdots p_k)^{-\ell} \mathbb{Q}[[\mathbf{p}]] $$

and a formal residue on it by the formula

$$ \mathrm{res}\left(\sum_{\mathbf{r} \in \mathbb{Z}^k} c_{\mathbf{r}} \mathbf{p}^{\mathbf{r}}\right) = c_{-1, \dots, -1}. \quad (21) $$

Lastly, we extend $h \mapsto \tilde{h}$ to a map from $\mathbb{Q}[\mathbf{p}]((t))$ to itself, $\mathcal{L}$ to a map from $\mathbb{Q}[\mathbf{p}]((t))$ to $\mathbb{Q}[\mathbf{p}^{-1}]((t))$, and $\mathrm{res}$ to a map from $\mathbb{Q}\langle\langle\mathbf{p}\rangle\rangle((t))$ to $\mathbb{Q}((t))$, by making each of those maps act coefficient-wise.

**Lemma 30.** For any polynomials $f$ and $g$ in $\mathbb{Q}[\mathbf{p}]$,

$$ \langle e^f, e^{tg} \rangle = \mathrm{res}(e^f \mathcal{L}(e^{t\tilde{g}})). $$

*Proof.* For $U \in \mathbb{Q}[[\mathbf{p}]]$ and for $\mathbf{r} \in \mathbb{N}^k$, $U\mathcal{L}(\tilde{\mathbf{p}}^{\mathbf{r}})$ is an element of $\mathbb{Q}\langle\langle\mathbf{p}\rangle\rangle$, so that using (20), (18), and (21) in order, we derive:

$$ \begin{aligned} \mathrm{res}(U\mathcal{L}(\tilde{\mathbf{p}}^{\mathbf{r}})) &= \mathrm{res}(U \times 1^{r_1} \cdots k^{r_k} r_1! p_1^{-r_1-1} \cdots r_k! p_k^{-r_k-1}) \\ &= z_{\mathbf{r}} \mathrm{res}(Up_1^{-r_1-1} \cdots p_k^{-r_k-1}) = z_{\mathbf{r}} \times [\mathbf{p}^{\mathbf{r}}]U = \langle U, \mathbf{p}^{\mathbf{r}} \rangle. \end{aligned} $$

This formula extends by linearity to $\langle U, h \rangle = \mathrm{res}(U\mathcal{L}(\tilde{h}))$ for any $h \in \mathbb{Q}[\mathbf{p}]$. Upon specializing to $U = e^f$ and $h = g^n/n!$ before taking series in $t$, this makes the informal integral formula provided in [23, end of 7.1] completely algebraic, in the form of the formula

$$ \langle e^f, e^{tg} \rangle = \sum_{\ell \ge 0} \langle e^f, g^\ell \rangle \frac{t^\ell}{\ell!} = \sum_{\ell \ge 0} \mathrm{res} \left( e^f \mathcal{L}(\tilde{g}^\ell) \right) \frac{t^\ell}{\ell!} = \mathrm{res} \left( e^f \mathcal{L}(e^{t\tilde{g}}) \right). \quad \square $$

The ring $\mathbb{Q}\langle\langle\mathbf{p}\rangle\rangle((t))$ is a $W_{\mathbf{p}}(t)$-module with the usual actions: $p_i$ acts by multiplication and $\partial_i$ by partial differentiation with respect to $p_i$. Let $K$ be the subspace of all elements of $\mathbb{Q}\langle\langle\mathbf{p}\rangle\rangle((t))$ that do not contain any monomial $\mathbf{p}^\mathbf{r}t^m$ with $r_1, \dots, r_k$ all negative. In other words,

$$ K = \sum_{i=1}^{k} \mathbb{Q}[[\mathbf{p}]] [p_1^{-1}, \dots, p_{i-1}^{-1}, p_{i+1}^{-1}, \dots, p_k^{-1}]((t)). $$

This subspace has the property to be a sub-$W_{\mathbf{p}}(t)$-module of $\mathbb{Q}\langle\langle\mathbf{p}\rangle\rangle((t))$ and to be contained in the kernel of the residue map $\mathbb{Q}\langle\langle\mathbf{p}\rangle\rangle((t)) \to \mathbb{Q}((t))$.

Now, let $f$ and $g$ be two polynomials in $\mathbb{Q}[\mathbf{p}]$. We provide in Lemma 31 an explicit construction of a derivation $\partial_t$ satisfying (19). We remark that this construction is simpler than the general approach presented in Section 4.3.
**Lemma 31.** The $W_x(t)$-linear map $L: a \mapsto a \tilde{g}(\frac{\partial f}{\partial p_1} - \partial_1, \dots, \frac{\partial f}{\partial p_k} - \partial_k)$ defines a derivation on $W_x(t)/S$ by

$$ \partial_t \cdot \mathrm{pr}_S(a) = \mathrm{pr}_S\left(\frac{\partial a}{\partial t} + L(a)\right). \qquad (22) $$

This derivation commutes with $p$ and $\partial_p$.

*Proof.* Let $\phi$ be the $W_p$-linear endomorphism of $W_p(t)$ defined by $\phi(a) = \frac{\partial a}{\partial t} + L(a)$. To show that $\partial_t$ is well-defined, it suffices to verify that $\phi(S) \subseteq S$. Consider the generators

$$ s_i := p_i - t \frac{\partial \tilde{g}}{\partial p_i} \left( \frac{\partial f}{\partial p_1} - \partial_1, \dots, \frac{\partial f}{\partial p_k} - \partial_k \right) $$

of $S$. Using the commutation rule $p_i(\frac{\partial f}{\partial p_j} - \partial_j) = (\frac{\partial f}{\partial p_j} - \partial_j)p_i + \delta_{i,j}$, one obtains for any polynomial $q(X_1, \dots, X_k)$

$$ \begin{aligned} p_i q\left(\frac{\partial f}{\partial p_1} - \partial_1, \dots, \frac{\partial f}{\partial p_k} - \partial_k\right) &= q\left(\frac{\partial f}{\partial p_1} - \partial_1, \dots, \frac{\partial f}{\partial p_k} - \partial_k\right) p_i + \frac{\partial q}{\partial X_i}\left(\frac{\partial f}{\partial p_1} - \partial_1, \dots, \frac{\partial f}{\partial p_k} - \partial_k\right). \end{aligned} $$

So, for each $i$, a computation shows that

$$ \phi(s_i) = \tilde{g}\left(\frac{\partial f}{\partial p_1} - \partial_1, \dots, \frac{\partial f}{\partial p_k} - \partial_k\right) s_i \in S. $$

Next, for any rational function $R(t)$ and any $i$, we check

$$ \phi(R(t)s_i) = R(t)\phi(s_i) + R'(t)s_i \in S, $$

by $\mathbb{Q}(t)$-linearity of $L$. Now, $S$ is generated as a $W_x$-module by the family $(R(t)s_i)_{R,i}$, so, by $W_x$-linearity of $L$, and thus of $\phi$, we get the inclusion $\phi(S) \subseteq S$. The $W_x(t)$-linearity of $L$ and the definition (22) imply, for any $R(t) \in \mathbb{Q}(t)$,

$$ \begin{aligned} \partial_t R(t) \cdot \mathrm{pr}_S(a) = \partial_t \cdot \mathrm{pr}_S(R(t)a)) = \mathrm{pr}_S\left(\frac{\partial(R(t)a)}{\partial t} + R(t)L(a)\right) = \\ R(t)\partial_t \cdot \mathrm{pr}_S(a) + R'(t) \mathrm{pr}_S(a). \end{aligned} $$

In other words, $\partial_t$ is a derivation. A similar but simpler calculation shows that it commutes with $p$ and $\partial_p$. $\square$

For the same polynomials $f$ and $g$, let $\Xi_{f,g}$ be the class of $e^f \mathcal{L}(e^{t\tilde{g}})$ modulo $K$.

**Lemma 32.** For any $a \in S$, the relation $a \cdot \Xi_{f,g} = 0$ holds. Moreover,

$$ \frac{\partial}{\partial t} \cdot \Xi_{f,g} = \tilde{g}\left(\frac{\partial f}{\partial p_1} - \partial_1, \dots, \frac{\partial f}{\partial p_k} - \partial_k\right) \cdot \Xi_{f,g}. $$
*Proof.* The definition (20) of the formal Laplace transform implies the following formulas, valid for $\mathbf{r} \in \mathbb{N}^k$:

$$ \mathcal{L}(p_i \cdot \mathbf{p}^r) = (r_i + 1)p_i^{-1} \mathcal{L}(\mathbf{p}^r) = -\partial_i \cdot \mathcal{L}(\mathbf{p}^r), \qquad (23) $$

$$ \mathcal{L}(\partial_i \cdot \mathbf{p}^r) = \mathcal{L}(r_i p_i^{-1} \mathbf{p}^r) = \begin{cases} p_i \cdot \mathcal{L}(\mathbf{p}^r), & \text{if } r_i \neq 0, \\ 0, & \text{otherwise.} \end{cases} \qquad (24) $$

In turn, for $h \in \mathbb{Q}[\mathbf{p}]$, this implies the formulas

$$ \mathcal{L}(p_i \cdot h) = -\partial_i \cdot \mathcal{L}(h), \qquad (25) $$

$$ \mathcal{L}(\partial_i \cdot h) = p_i \cdot \mathcal{L}(h) - p_i \cdot \mathcal{L}(h|_{p_i=0}). \qquad (26) $$

The last formula is not convenient because of the term involving $h|_{p_i=0}$. Fortunately, this term is in $K$, so we have the nicer formula

$$ \mathcal{L}(\partial_i \cdot h) \equiv p_i \cdot \mathcal{L}(h) \pmod{K}. \qquad (27) $$

Moreover, we have for any $h \in \mathbb{Q}\langle\langle\mathbf{p}\rangle\rangle((t))$,

$$ (\frac{\partial f}{\partial p_i} - \partial_i) \cdot e^f h = -e^f \partial_i \cdot h \qquad (28) $$

Therefore, we have

$$ \begin{aligned} & \left( p_i - t \frac{\partial \tilde{g}}{\partial p_i} \left( \frac{\partial f}{\partial p_1} - \partial_1, \dots, \frac{\partial f}{\partial p_k} - \partial_k \right) \right) \cdot e^f \mathcal{L}(e^{t\tilde{g}}) \qquad (29) \\ & = e^f \left( p_i - t \frac{\partial \tilde{g}}{\partial p_i} (-\partial_1, \dots, -\partial_k) \right) \cdot \mathcal{L}(e^{t\tilde{g}}), \quad \text{using (28)} \\ & \equiv e^f \mathcal{L} \left( \left( \partial_i - t \frac{\partial \tilde{g}}{\partial p_i} \right) \cdot e^{t\tilde{g}} \right) \quad \text{using (25) and (27)} \\ & \equiv e^f \mathcal{L}(0) \equiv 0 \pmod{K}. \end{aligned} $$

This proves the first statement about all $a \in S$ by $W_p(t)$-linearity. The second statement is proved similarly, starting with $\frac{\partial}{\partial t} - \tilde{g}(\frac{\partial f}{\partial p_1} - \partial_1, \dots, \frac{\partial f}{\partial p_k} - \partial_k)$ instead of the operator in (29). $\square$

**Lemma 33.** The $W_p(t)$-module $W_p(t)/S$ is holonomic.

*Proof.* Let $\tau$ be the automorphism of the $\mathbb{Q}(t)$-algebra $W_p(t)$ defined by

$$ \tau(p_i) = \frac{\partial f}{\partial p_i} - \partial_i \quad \text{and} \quad \tau(\partial_i) = p_i - t\tau\left(\frac{\partial \tilde{g}}{\partial p_i}\right). $$

(Note that the definition is not recursive since $\frac{\partial \tilde{g}}{\partial p_i}$ is a polynomial in $\mathbf{p}$. Also, $p_i$ commutes with $\frac{\partial \tilde{g}}{\partial p_i}$, making their images under $\tau$ commute as well. This justifies $\tau(\partial_i)\tau(p_i) = \tau(p_i)\tau(\partial_i) + 1$ to have an algebra morphism.) The inverse morphism is given by

$$ \tau^{-1}(p_i) = \partial_i + t \frac{\partial \tilde{g}}{\partial p_i} \quad \text{and} \quad \tau^{-1}(\partial_i) = \tau^{-1}\left(\frac{\partial f}{\partial p_i}\right) - p_i. $$
By definition, S is generated by $\tau(\partial_1), \dots, \tau(\partial_k)$. Therefore, as $W_p(t)$-modules,

$$W_p(t)/S \simeq W_p(t)/(W_p(t)\partial_1 + \dots + W_p(t)\partial_k) \simeq \mathbb{Q}(t)[\mathbf{p}].$$

Since the latter is holonomic, so is $W_p(t)/S$. $\square$

Let $M$ denote the holonomic module $W_p(t)/S$. All in all, we have a commutative diagram of $\mathbb{Q}(t)$-linear spaces, with all arrows commuting with the derivation $\partial_t$

$$\begin{tikzcd}
M \arrow[r] \arrow[d] & W_p(t) \cdot \Xi_{f,g} \arrow[d, "\text{res}"] \arrow[r, "\mathbb{Q}(\langle\langle \mathbf{p} \rangle\rangle((t)) / K"] & \\
\frac{M}{\partial M} \arrow[r] & W_p(t) \cdot \langle e^f, e^{tg} \rangle \arrow[r, "\text{res}"] & \mathbb{Q}((t)).
\end{tikzcd}$$

The class of 1 in $M$ is mapped to $\langle e^f, e^{tg} \rangle$ in $\mathbb{Q}((t))$. What Algorithm 4 computes, is an operator $L(t, \partial_t)$ such that $L(t, \partial_t) \cdot 1 \in \partial M$. This implies that $L(t, \partial_t) \cdot \langle e^f, e^{tg} \rangle = 0$.

*Remark 34.* At this point we can make explicit our remark that earlier creative-telescoping algorithms could not deal with our integrals. From the explicit definition of $g$ in the formula $R_k(t) = \langle e^f, e^{tg} \rangle$, we can prove that $\tilde{g}$ is always in the form $p_k + h(p_1, \dots, p_{k-1})$, for some polynomial $h$. So $p_k - t$ is always in $S$. Any integration algorithm that would work over $\mathbb{Q}(t, p_k)$, as many that are designed to apply to functions, would therefore consider that 1 is in the annihilator of the function to be integrated, which would lead to a wrong result.

## 6.2 Experimental results

We consider graph models that are either some model of $k$-regular (simple) graphs or some generalization with loops and/or multiple edges and/or degrees in the set $\{1, 2, \dots, k\}$ instead of $\{k\}$. Given such a graph model, the theory in [22] provides immediate formulas for the polynomials $f$ and $g$. Obtaining the ideal of the lemma is easily implemented as a simple non-commutative substitution. To this end, we used Maple's OreAlgebra package. After converting<fnref n="2" /> from Maple notation to the notation of MultivariateCreativeTelescoping.jl, we could use the latter to obtain the wanted ODEs, appealing to the implementation of our optimized Algorithm 5. This was done for $2 \le k \le 8$ and the degree sets $K = \{k\}$ and $K = \{1, 2, \dots, k\}$. We collected the computed ODEs and made them available on the web<fnref n="3" />. For $k \le 7$, they are the same as with the Maple implementation that accompanies Chyzak and Mishna's article [22].

<footnote n="2">²Maple uses its commutative product `*` to represent monomials, so that both Maple inputs `t*dt` and `dt*t` represent the element $t\partial_t$. This is no problem inside Maple, where non-commutative products are computed by the command `skew_product`. But naively serializing an operator from Maple by `lprint` can lead to strings with a different interpretation in Brochet's Julia implementation. We automated a rewrite of those strings to move all derivatives to the right.</footnote>

<footnote n="3">³See https://files.inria.fr/chyzak/kregs/.</footnote>
<table><thead><tr><th colspan="3">graph</th><th colspan="2">ODE</th><th colspan="4">time and memory</th></tr><tr><th>k</th><th>l</th><th>e</th><th>ord</th><th>deg</th><th>f5</th><th>mct</th><th>total</th><th>rss</th></tr></thead><tbody><tr><td>2</td><td>ll</td><td>se</td><td>1</td><td>2</td><td>0.04</td><td>0.05</td><td>18</td><td>0.63</td></tr><tr><td>3</td><td>ll</td><td>se</td><td>2</td><td>11</td><td>0.04</td><td>0.05</td><td>18</td><td>0.62</td></tr><tr><td>4</td><td>ll</td><td>se</td><td>2</td><td>14</td><td>0.05</td><td>0.05</td><td>19</td><td>0.62</td></tr><tr><td>5</td><td>ll</td><td>se</td><td>6</td><td>125</td><td>0.05</td><td>0.59</td><td>17</td><td>0.65</td></tr><tr><td>6</td><td>ll</td><td>se</td><td>6</td><td>145</td><td>0.23</td><td>1.0</td><td>20</td><td>0.66</td></tr><tr><td>7</td><td>ll</td><td>se</td><td>20</td><td>1683</td><td>8.6</td><td>303</td><td>330</td><td>4.6</td></tr><tr><td>7</td><td>ll</td><td>me</td><td>20</td><td>1683</td><td>8.4</td><td>300</td><td>326</td><td>4.3</td></tr><tr><td>7</td><td>la</td><td>se</td><td>20</td><td>1683</td><td>8.3</td><td>301</td><td>328</td><td>4.3</td></tr><tr><td>7</td><td>la</td><td>me</td><td>20</td><td>1683</td><td>8.3</td><td>299</td><td>326</td><td>4.4</td></tr><tr><td>7</td><td>lh</td><td>se</td><td>20</td><td>1683</td><td>8.6</td><td>310</td><td>337</td><td>5.9</td></tr><tr><td>7</td><td>lh</td><td>me</td><td>20</td><td>1683</td><td>8.3</td><td>322</td><td>349</td><td>5.8</td></tr><tr><td>8</td><td>ll</td><td>se</td><td>19</td><td>1793</td><td>244</td><td>832</td><td>1095</td><td>6.5</td></tr><tr><td>8</td><td>ll</td><td>me</td><td>19</td><td>1793</td><td>247</td><td>831</td><td>1097</td><td>6.7</td></tr><tr><td>8</td><td>la</td><td>se</td><td>19</td><td>1793</td><td>244</td><td>831</td><td>1094</td><td>6.0</td></tr><tr><td>8</td><td>la</td><td>me</td><td>19</td><td>1793</td><td>244</td><td>829</td><td>1093</td><td>6.0</td></tr><tr><td>8</td><td>lh</td><td>se</td><td>35</td><td>6204</td><td>393</td><td>23069</td><td>23481</td><td>5.4</td></tr><tr><td>8</td><td>lh</td><td>me</td><td>35</td><td>6200</td><td>393</td><td>23111</td><td>23524</td><td>5.4</td></tr></tbody></table>

<table><caption>Table 1: Computation of ODEs for k-regular graph models $K = \{k\}$ (left) and $K = \{1, \dots, k\}$ (right). A graph model is input by the triple $(k, l, e)$, where $l$ and $e$ describe if loops and edges are allowed. For each output ODE, the order ('ord') and the coefficient degree ('deg') are given. All times are given in seconds. The two main computation steps are preparing a Gröbner basis ('f5') and running Algorithm 5 on it ('mct'). The maximum memory used is given in GB ('rss'). See details in Section 6.2.</caption><thead><tr><th colspan="3">graph</th><th colspan="2">ODE</th><th colspan="4">time and memory</th></tr><tr><th>k</th><th>l</th><th>e</th><th>ord</th><th>deg</th><th>f5</th><th>mct</th><th>total</th><th>rss</th></tr></thead><tbody><tr><td>2</td><td>ll</td><td>se</td><td>1</td><td>3</td><td>0.04</td><td>0.04</td><td>17</td><td>0.62</td></tr><tr><td>3</td><td>ll</td><td>se</td><td>2</td><td>11</td><td>0.04</td><td>0.13</td><td>18</td><td>0.64</td></tr><tr><td>4</td><td>ll</td><td>se</td><td>3</td><td>29</td><td>0.04</td><td>0.07</td><td>18</td><td>0.62</td></tr><tr><td>5</td><td>ll</td><td>se</td><td>6</td><td>125</td><td>0.06</td><td>0.69</td><td>19</td><td>0.65</td></tr><tr><td>6</td><td>ll</td><td>se</td><td>10</td><td>425</td><td>0.31</td><td>11</td><td>31</td><td>0.86</td></tr><tr><td>7</td><td>ll</td><td>se</td><td>20</td><td>1683</td><td>8.3</td><td>316</td><td>343</td><td>5.8</td></tr><tr><td>7</td><td>ll</td><td>me</td><td>20</td><td>1683</td><td>8.5</td><td>321</td><td>348</td><td>5.8</td></tr><tr><td>7</td><td>la</td><td>se</td><td>20</td><td>1683</td><td>8.7</td><td>324</td><td>352</td><td>5.8</td></tr><tr><td>7</td><td>la</td><td>me</td><td>20</td><td>1683</td><td>8.8</td><td>322</td><td>349</td><td>5.8</td></tr><tr><td>7</td><td>lh</td><td>se</td><td>20</td><td>1683</td><td>8.7</td><td>310</td><td>337</td><td>5.6</td></tr><tr><td>7</td><td>lh</td><td>me</td><td>20</td><td>1683</td><td>8.2</td><td>323</td><td>349</td><td>5.8</td></tr><tr><td>8</td><td>ll</td><td>se</td><td>35</td><td>6201</td><td>389</td><td>23198</td><td>23605</td><td>5.4</td></tr><tr><td>8</td><td>ll</td><td>me</td><td>35</td><td>6200</td><td>386</td><td>23586</td><td>23991</td><td>5.4</td></tr><tr><td>8</td><td>la</td><td>se</td><td>35</td><td>6204</td><td>401</td><td>23495</td><td>23915</td><td>5.4</td></tr><tr><td>8</td><td>la</td><td>me</td><td>35</td><td>6205</td><td>393</td><td>23188</td><td>23600</td><td>5.4</td></tr><tr><td>8</td><td>lh</td><td>se</td><td>35</td><td>6205</td><td>387</td><td>22745</td><td>23152</td><td>5.5</td></tr><tr><td>8</td><td>lh</td><td>me</td><td>35</td><td>6205</td><td>394</td><td>23440</td><td>23853</td><td>5.4</td></tr></tbody></table>

To the best of our knowledge, the ODEs for $k = 8$ are obtained for the first time. Table 1 displays some parameters related to the calculations and to the results we obtained. There, each graph model is described by a triple $(k, l, e)$ where: $k$ determines the set of allowed degrees, either by $K = \{k\}$ for the left part of Table 1 or by $K = \{1, \dots, k\}$ for the right part; $l$ is one of 'll' for loop-less graphs, 'la' for graphs with loops allowed and contributing 2 to the degree, and 'lh' for graphs with loops allowed and contributing 1 to the degree; $e$ is either 'se' for graphs with simple edges or 'me' for graphs with multiple edges allowed. For each graph model, the resulting ODE has order and degree given in columns 'ord' and 'deg'. The total time for computing the ODE is given in column 'total' and decomposes as follows: the time to prepare generators for the ideal $S$ is negligible; the time to make a Gröbner basis out of them is given in column 'f5'; the time to extract next the ODE by our new algorithm is given in column 'mct'. We note that the time for 'mct' always dominates. The total time includes the compilation time, which explains why its value exceeds the sum of 'f5' and 'mct'. The maximal peak of memory usage is listed in column 'rss'.

By way of comparison, the model (se, ll, 7) could be computed by the method
and implementation of [22] in 3.22 · 10⁴ seconds (almost 9 hours), which is
roughly 100 times as much as the 330 seconds (5.5 minutes) needed by our
new algorithm and implementation. This can partly be explained by the lack
of efficient evaluation/interpolation methods in the implementation of [22] and
by the choice in [22] to continue the calculation by factoring the polynomial
coefficients of the ODE.

**Acknowledgements.** This work has been supported by the European Research Council under the European Union's Horizon Europe research and innovation programme, grant agreement 101040794 (10000 DIGITS).

References

[1] W. W. Adams and P. Loustaunau. *An Introduction to Gröbner Bases*. American Mathematical Society, July 1994.

[2] R. Bahloul. "Gröbner bases in D-modules: application to Bernstein-Sato polynomials". In: *Two algebraic byways from differential equations: Gröbner bases and quivers*. Vol. 28. Algorithms Comput. Math. Springer, Cham, 2020, pp. 75-93.

[3] T. Becker and V. Weispfenning. *Gröbner Bases*. Springer New York, 1993.
DOI: 10/cfhwn9.

[4] I. N. Bernshtein. "Modules over a ring of differential operators. Study of the fundamental solutions of equations with constant coefficients - Functional Analysis and Its Applications". In: (1971). DOI: 10/fcjn47.

[5] F. Beukers. "Irrationality of π², periods of an elliptic curve and Γ₁(5)". In: 31. Progr. Math. Mass. (1983), pp. 47–66.

[6] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah. "Julia: A Fresh Approach to Numerical Computing". In: *SIAM Rev.* 59.1 (2017), pp. 65–98. DOI: 10/f9wkpj.

[7] J.-E. Björk. *Rings of differential operators*. Vol. 21. North-Holland Mathematical Library. North-Holland Publishing Co., Amsterdam-New York, 1979, pp. xvii+374.

[8] A. Borel, P.-P. Grivel, B. Kaup, A. Haefliger, B. Malgrange, and F. Ehlers. *Algebraic D-modules*. Vol. 2. Perspectives in Mathematics. Academic Press, Inc., Boston, MA, 1987, pp. xii+355.

[9] A. Bostan, S. Chen, F. Chyzak, and Z. Li. "Complexity of Creative Tele-
scoping for Bivariate Rational Functions". In: *Proc. 35th Int. Symp. Symb.
Algebr. Comput*. ACM, 2010, pp. 203–210. DOI: 10/d4tf3r.

[10] A. Bostan, S. Chen, F. Chyzak, Z. Li, and G. Xin. "Hermite Reduction and Creative Telescoping for Hyperexponential Functions". In: *Proc. 38th Int. Symp. Symb. Algebr. Comput*. ISSAC '13. ACM, 2013, pp. 77–84.
DOI: 10/ggck9v.
[11] A. Bostan, F. Chyzak, P. Lairez, and B. Salvy. "Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions". In: *ISSAC'18*. Ed. by É. Schost. ACM Press, 2018, pp. 95–102. DOI: 10/ddv8.

[12] A. Bostan, P. Lairez, and B. Salvy. "Creative Telescoping for Rational Functions Using the Griffiths–Dwork Method". In: *Proc. ISSAC 2013*. ACM, 2013, pp. 93–100. DOI: 10/ggcmbk.

[13] F. Castro. "Calculs effectifs pour les idéaux d'opérateurs différentiels". In: *Géométrie algébrique et applications, III (La Rábida, 1984)*. Vol. 24. Travaux en Cours. Hermann, Paris, 1987, pp. 1–19. DOI: 10/d2zr4n.

[14] F. Castro. "Théorème de division pour les opérateurs différentiels et calcul des multiplicités". thèse de 3e cycle. Université Paris VII, 1984.

[15] S. Chen, L. Du, and M. Kauers. "Hermite Reduction for D-finite Functions via Integral Bases". In: *Proc. 2023 Int. Symp. Symb. Algebr. Comput. ISSAC '23*. Association for Computing Machinery, July 24, 2023, pp. 155–163. DOI: 10/g865z3.

[16] S. Chen, M. van Hoeij, M. Kauers, and C. Koutschan. "Reduction-Based Creative Telescoping for Fuchsian D-finite Functions". In: *J. Symb. Comput*. 85 (Mar. 1, 2018), pp. 108–127. DOI: 10/ggck9k.

[17] S. Chen, H. Huang, M. Kauers, and Z. Li. "A Modified Abramov-Petkovšek Reduction and Creative Telescoping for Hypergeometric Terms". In: *Proc. 40th Int. Symp. Symb. Algebr. Comput*. Ed. by ACM. 2015. DOI: 10/f3s5fq.

[18] S. Chen, M. Kauers, and C. Koutschan. "Reduction-Based Creative Telescoping for Algebraic Functions". In: *Proc. ACM Int. Symp. Symb. Algebr. Comput*. ISSAC '16. ACM, 2016, pp. 175–182. DOI: 10/ggck9z.

[19] F. Chyzak and B. Salvy. "Non-commutative Elimination in Ore Algebras Proves Multivariate Identities". In: *Journal of Symbolic Computation* 26.2 (1998), pp. 187–227. DOI: 10/bzg7x8.

[20] F. Chyzak. "An Extension of Zeilberger's Fast Algorithm to General Holo- nomic Functions". In: *Discrete Math*. 217.1–3 (2000), pp. 115–134. DOI: 10/drkkn6.

[21] F. Chyzak. "The ABC of Creative Telescoping: Algorithms, Bounds, Com- plexity". 64 pages. Memoir of accreditation to supervise research (HDR). Université d'Orsay, Apr. 2014.

[22] F. Chyzak and M. Mishna. "Differential equations satisfied by generat- ing functions of 5-, 6-, and 7-regular labelled graphs: a reduction-based approach". Preprint. June 2024.

[23] F. Chyzak, M. Mishna, and B. Salvy. "Effective scalar products of D- finite symmetric functions". In: *Journal of Combinatorial Theory. Series A* 112.1 (2005), pp. 1–43.
[24] S. C. Coutinho. *A Primer of Algebraic D-Modules*. London Mathematical Society Student Texts. Cambridge University Press, 1995.

[25] D. Cox, J. Little, and D. O'Shea. *Using algebraic geometry*. Vol. 185. Graduate Texts in Mathematics. Springer-Verlag, New York, 1998, pp. xii+499.
DOI: 10/ghmnj2.

[26] A. Dimca. "On the Milnor Fibrations of Weighted Homogeneous Polynomials". In: *Compositio Math.* 76.1-2 (1990), pp. 19-47.

[27] B. Dwork. "On the zeta function of a hypersurface". en. In: *Publications Mathématiques de l'IHÉS* 12 (1962), pp. 5-68.

[28] B. Dwork. "On the Zeta Function of a Hypersurface: II". In: *Annals of Mathematics* 80.2 (1964), pp. 227-299.

[29] J.-C. Faugère. "A new efficient algorithm for computing Gröbner bases (F4)". In: *Journal of Pure and Applied Algebra* 139.1 (1999), pp. 61-88.
DOI: 10/bpq5dx.

[30] A. Galligo. "Some algorithmic questions on ideals of differential operators". In: *Eurocal'85, Vol. 2 (Linz, 1985)*. Vol. 204. Lecture Notes in Comput. Sci. Berlin: Springer, 1985, pp. 413-421.

[31] J. von zur Gathen and J. Gerhard. *Modern Computer Algebra*. Cambridge University Press, 1999.

[32] P. A. Griffiths. "On the Periods of Certain Rational Integrals: I". In: *Annals of Mathematics* 90.3 (1969), pp. 460-495.
DOI: 10/dq3zc6.

[33] J. van der Hoeven. "Constructing Reductions for Creative Telescoping". In: *AAECC* 32.5 (Nov. 1, 2021), pp. 575-602.
DOI: 10/g865z4.

[34] A. Kandri-Rody and V. Weispfenning. "Noncommutative Gröbner bases in algebras of solvable type". In: *Journal of Symbolic Computation* 9.1 (1990), pp. 1-26.

[35] M. Kauers, M. Jaroschek, and F. Johansson. "Ore polynomials in Sage". In: *Computer Algebra and Polynomials: Applications of Algebra and Number Theory*. Ed. by J. Gutierrez, J. Schicho, and M. Weimann. Cham: Springer International Publishing, 2015, pp. 105-125.
DOI: 10/f3s5ht.

[36] C. Koutschan. "A Fast Approach to Creative Telescoping". In: *Math. Comput. Sci.* 4.2-3 (2010), pp. 259-266.
DOI: 10/bhb6sv.

[37] C. Koutschan. *Holonomic Functions, User's Guide*. 10-01. RISC Report Series, University of Linz, Austria, 2010.

[38] P. Lairez. "Computing Periods of Rational Integrals". In: *Math. Comput.* 85.300 (2016), pp. 1719-1752.
DOI: 10/ggck95.

[39] P. Lairez. "Computing periods of rational integrals". In: *Mathematics of Computation* 85.300 (2016), pp. 1719-1752.
DOI: 10/ggck95.

[40] P. Lairez. "Axioms for a theory of signature bases". In: *Journal of Symbolic Computation* 123 (2024), p. 102275.
DOI: 10/k99k.
[41] S. Laporta. "High-Precision Calculation of Multiloop Feynman Integrals by Difference Equations". In: *Int. J. Mod. Phys. A* 15.32 (Dec. 30, 2000), pp. 5087-5159. DOI: 10/czhxfg.

[42] V. Levandovskyy and H. Schönenmann. "Plural: a computer algebra system for noncommutative polynomial algebras". In: *Proceedings of the 2003 international symposium on Symbolic and algebraic computation - ISSAC'03*. ACM Press, 2003. DOI: 10/b593w7.

[43] T. Oaku. "Algorithms for Integrals of Holonomic Functions over Domains Defined by Polynomial Inequalities". In: *J. Symb. Comput.* 50 (2013), pp. 1-27. DOI: 10/f4hpj4.

[44] T. Oaku and N. Takayama. "Algorithms for D-modules — Restriction, Tensor Product, Localization, and Local Cohomology Groups". In: *Journal of Pure and Applied Algebra* 156.2 (Feb. 2001), pp. 267-308. DOI: 10/bct97n.

[45] T. Oaku, N. Takayama, and U. Walther. "A localization algorithm for D-modules". In: *Journal of Symbolic Computation* 29.4-5 (2000). Special issue on Symbolic Computation in Algebra, Analysis, and Geometry (Berkeley, CA, 1998), pp. 721-728. DOI: 10/fgnnmm.

[46] M. Saito, B. Sturmfels, and N. Takayama. *Gröbner deformations of hypergeometric differential equations*. Vol. 6. Algorithms and Computation in Mathematics. Berlin: Springer-Verlag, 2000, pp. viii+254.

[47] Y. Sun, D. Wang, X. Ma, and Y. Zhang. "A signature-based algorithm for computing Gröbner bases in solvable polynomial algebras". In: *Proceedings of the 37th International Symposium on Symbolic and Algebraic Computation*. ISSAC '12. Grenoble, France: Association for Computing Machinery, 2012, pp. 351-358. DOI: 10/ghtkmx.

[48] N. Takayama. "Gröbner basis and the problem of contiguous relations". In: *Japan Journal of Applied Mathematics* 6.1 (1989), pp. 147-160. DOI: 10/fqjhjg.

[49] N. Takayama. "An algorithm of constructing the integral of a module: an infinite dimensional analog of Gröbner basis". In: *Symbolic and Algebraic Computation*. Proceedings of ISSAC'90 (Kyoto, Japan). ACM and Addison-Wesley, 1990, pp. 206-211.

[50] N. Takayama. "Gröbner Basis, Integration and Transcendental Functions". In: *Proc. ISSAC 1990* (Tokyo, Japan). ACM, 1990, pp. 152-156. DOI: 10/d62qqj.

[51] C. Traverso. "Gröbner trace algorithms". In: *Symbolic and Algebraic Computation*. Ed. by P. Gianni. Berlin, Heidelberg: Springer Berlin Heidelberg, 1989, pp. 125-138.

[52] H. S. Wilf and D. Zeilberger. "An algorithmic proof theory for hypergeometric (ordinary and "q") multisum/integral identities". In: *Invent. Math.* 108.3 (1992), pp. 575-633. DOI: 10/fnvh8h.
[53] D. Zeilberger. "A Fast Algorithm for Proving Terminating Hypergeometric Identities". In: *Discrete Mathematics* 80.2 (Mar. 1, 1990), pp. 207–211.
DOI: 10/crrzhp