# Minimization of differential equations
and
algebraic values of $E$-functions

A. Bostan  Inria, Université Paris-Saclay, 1 rue Honoré d’Estienne d’Orves,
91120 Palaiseau, France  ,  T. Rivoal  Institut Fourier, CNRS et Université
Grenoble Alpes, CS 40700, 38058 Grenoble cedex 9  and  B. Salvy  Univ Lyon,
EnsL, UCBL, CNRS, Inria, LIP, F-69342, LYON Cedex 07, France

###### Abstract.

A power series being given as the solution of a linear differential equation
with appropriate initial conditions, minimization consists in finding a non-
trivial linear differential equation of minimal order having this power series
as a solution. This problem exists in both homogeneous and inhomogeneous
variants; it is distinct from, but related to, the classical problem of
factorization of differential operators. Recently, minimization has found
applications in Transcendental Number Theory, more specifically in the
computation of non-zero algebraic points where Siegel’s $E$-functions take
algebraic values. We present algorithms for these questions and discuss
implementation and experiments.

###### Key words and phrases:

Linear differential operators, Minimization, Factorization, Desingularization,
$E$-functions, Siegel-Shidlovskii theorem, Beukers’ algorithm, Adamczewski-
Rivoal algorithm

###### 2010 Mathematics Subject Classification:

68W30, 11J81, 16S32, 34M15, 33F10

##  1\. Introduction

###  1.1. Minimization

A linear differential equation

${\mathcal{L}}(y(z))=a_{r}(z)y^{(r)}(z)+\dots+a_{0}(z)y(z)=0$ (1)

with polynomial coefficients in ${\mathbb{Q}}[z]$ is given, together with
initial conditions specifying a unique formal power series solution
$S\in{\mathbb{Q}}[[z]]$. In its homogeneous variant, the problem of
minimization is to find a homogeneous linear differential equation of minimal
order and with polynomial coefficients in ${\mathbb{Q}}[z]$ having $S$ as a
solution. In the inhomogeneous version, the input is the same, but in the
output, a non-zero polynomial right-hand side in ${\mathbb{Q}}[z]$ is also
possible, which may allow for the existence of an equation of even smaller
order. Both these problems exist for other fields of coefficients and a large
part of our discussion extends to such situations; we focus here on the case
of the field ${\mathbb{Q}}$ to keep the discussion simple and reflect more
closely the capabilities of our implementation.

When the origin is an _ordinary point_ , i.e. $a_{r}(0)\neq 0$, initial
conditions are given as the values of $(S(0),\dots,S^{(r-1)}(0))$. Otherwise,
the origin is a _singularity_ of the equation (1). It may still have power
series solutions. The definition of initial conditions in this case is more
delicate. It is discussed in Section 2.1.

###  1.2. Relation to factorization of linear differential operators

The problems of factorization and minimization are closely related since they
are both concerned with finding (right) factors of linear differential
operators. But they are different problems. For instance, an equation can be
minimal even when the operator factors. A simple example is given by the
equation $(1-z)y^{\prime\prime}-y^{\prime}=0$ and its power series solution
$S=\ln(1-z)$. The corresponding operator clearly has
$\partial_{z}:=\frac{d}{dz}$ as a right factor of order 1, but no homogeneous
equation of order 1 can have a solution with a logarithmic singularity. Also,
in general, a linear differential operator has infinitely many factorizations
and the problem of minimization is to find a minimal, not necessarily
irreducible, right factor that vanishes at the solution $S$. Thus one cannot
simply use an existing implementation of a factorization algorithm in order to
solve the minimization problem.

Still, factorization and minimization share many algorithmic tools. Indeed,
the algorithm we present in Section 2 is obtained by combining sub-algorithms
of van Hoeij’s description of his factorization algorithm [50, 51], exploiting
the fact that the situation of minimization is made easier by the extra
information provided by the input power series $S$. One of our contributions
lies in a pedagogical presentation of the necessary tools, unencumbered by the
intricacies of the general factorization problem. A simple example of the
difference of behaviour between the minimization problem and the general
factorization problem is illustrated in Section 4.1. Furthermore, in
applications, it is sometimes possible to take advantage of further structure
that the minimal operator is known to possess, such as being Fuchsian, or
having at most one irregular singularity, at infinity, as happens for
$E$-functions (see Section 3).

###  1.3. Applications to $E$-functions

The above algorithmic considerations have immediate applications in number
theory, more precisely to the transcendence theory of the values taken at
algebraic points by $E$-functions. These power series in
${\overline{\mathbb{Q}}}[[z]]$ are solutions of linear differential equations
over ${\overline{\mathbb{Q}}}(z)$; they have been introduced and studied by
Siegel in 1929 as a generalization of the exponential function. We refer to
Section 3 for their definition, statements of results and bibliographic
references. $E$-functions have been studied in depth by Siegel, Shidlovskii,
Nesterenko, André, Beukers and others. Adamczewski and Rivoal have given an
algorithm [5] that determines the finite list of algebraic numbers $\alpha$
such that $f(\alpha)$ is also algebraic, given as input an $E$-function $f(z)$
(represented by a differential equation and initial conditions). The first two
steps of this algorithm rely on the computation of a minimal homogeneous
equation and of a minimal inhomogeneous equation for $f$. The implementation
of that algorithm is not efficient in practice. The new version presented here
enables to speed up the necessary computations. The final step is based on a
process of desingularization of differential systems, due to Beukers, and of
independent interest, that we make more explicit.

Similar questions have been addressed for the class of Mahler functions over
${\overline{\mathbb{Q}}}(z)$, i.e., formal power series in
${\overline{\mathbb{Q}}}[[z]]$ solutions of an equation
$\sum_{j=0}^{d}p_{j}(z)f(z^{r^{j}})=0$ where
$p_{j}(z)\in{\overline{\mathbb{Q}}}[z]$. An analogue of Beukers’ lifting
theorem has been proved by Adamczewski and Faverjon [4] and by Philippon [38].
As a consequence, there exists an algorithm that determines the list of
algebraic numbers $\alpha$ such that $f(\alpha)$ is algebraic, for any given
Mahler function $f(z)$. It seems however that the effective implementation of
the underlying algorithm has not yet been done.

### Structure of the article

In Section 2 we describe our new minimization algorithm, both in its
homogeneous form (Section 2.2) and in its inhomogeneous variant (Section 2.4).
Both crucially rely on computations of degree bounds which are described in
Section 2.3. In Section 3 we discuss an application of the minimization
algorithm to the problem of computing the set of algebraic values taken by
$E$-functions over ${\mathbb{Q}}$ at algebraic points. The new algorithm is a
practical variant of the Adamczewski-Rivoal algorithm recalled in Section 3.1,
itself based on Beukers’ desingularization procedure described and enhanced in
Section 3.2. The algorithm from Section 3.1 is then applied in Section 3.3 to
obtain an effective decomposition of an $E$-function over ${\mathbb{Q}}$ as a
polynomial plus a polynomial multiple of a purely transcendental $E$-function.
Extensions to $E$-functions with coefficients in a number field, respectively
to $E$-functions in Siegel’s original sense are discussed in Sections 3.4 and
3.5. Finally, Section 4 describes our implementation of the algorithms and
illustrates it with a few examples and timings. In particular, two infinite
families of $E$-functions that take algebraic values at non-trivial algebraic
points are presented in Section 4.3.

##  2\. Minimization algorithm

###  2.1. Power series solutions

We recall properties of linear differential equations that can be found in the
classical treatises of Ince [28, Chap. XVI, XVII] or Poole [39, Chap. V].
Moreover, the presentation is specialized to the case of coefficients $a_{i}$
of Eq. 1 that are polynomials rather than formal power series.

The image by $\mathcal{L}$ of a monomial $z^{s}$ with $s\in\mathbb{N}$ is a
polynomial

$f(s,z)=z^{s+g}(p_{0}(s)+p_{1}(s)z+\dots+p_{t}(s)z^{t}),\qquad-r\leq g,\quad
0\leq t,$ (2)

with polynomials $p_{i}(s)$ of degree at most $r$ whose coefficients depend on
those of the $a_{i}$ and $p_{0}\neq 0$. The polynomial $p_{0}$ is called the
_indicial polynomial_ of $\mathcal{L}$ at 0. By linear combination, the image
by $\mathcal{L}$ of a formal power series $S(z)=\sum_{i\geq 0}c_{i}z^{i}$ is
the formal power series

$\mathcal{L}(S)=\sum_{i\geq 0}c_{i}f(i,z).$

The coefficients of $z^{k}$ for $k=g,g+1,\dots$ in $\mathcal{L}(S)=0$ give the
equations

$c_{0}p_{0}(0)=0,\ c_{0}p_{1}(0)+c_{1}p_{0}(1)=0,\dots,\
c_{0}p_{t-1}(0)+\dots+c_{t-1}p_{0}(t-1)=0,$ (3)

and the linear recurrence of order $t$

$c_{i}p_{t}(i)+\dots+c_{t+i}p_{0}(t+i)=0,\quad i\geq 0.$ (4)

These equations imply that the valuation of $S$ (the index of its first non-
zero coefficient) is a zero of the indicial polynomial $p_{0}$. Let

$\mathcal{Z}_{\mathcal{L}}=\\{k\in\mathbb{N}\mid p_{0}(k)=0\\}$

be the set of nonnegative integer roots of the indicial polynomial of
$\mathcal{L}$ at 0. For all $i\not\in\mathcal{Z}_{\mathcal{L}}$, the
coefficient $c_{i}$ is determined from the previous ones by the $(i+1)$th
equation of the infinite system (3)–(4). For this reason, the _initial
conditions_ of the differential equation (1) are the values of $y^{(i)}(0)$
for $i\in\mathcal{Z}_{\mathcal{L}}$, as all the other ones are determined by
the system (3)–(4). In the non-singular case when $a_{r}(0)\neq 0$, the
indicial polynomial is $p_{0}(s)=s(s-1)\dotsm(s-r+1)$ and then
$\mathcal{Z}_{\mathcal{L}}=\\{0,1,\dots,r-1\\}$, recovering the usual
definition. This discussion leads to the following result that will be used to
find right factors of $\mathcal{L}$. (See Prop. 4.3 and Section 4.3 in [16]
for similar considerations.)

######  Lemma 2.1.

With the notation above, let $S$ be a power series solution of $\mathcal{L}$
and $\mathcal{M}$ be a right factor of $\mathcal{L}$. If there exists a
polynomial $T$ such that $T^{(i)}(0)=S^{(i)}(0)$ for all
$i\in\mathcal{Z}_{\mathcal{L}}$ and
$\mathcal{M}(T)=O(z^{\max\mathcal{Z}_{\mathcal{L}}+1})$, then
$\mathcal{M}(S)=0$.

###### Proof.

Applying the discussion above to $\mathcal{M}$ shows that the coefficients of
$T$ satisfy the first $\max\mathcal{Z}_{\mathcal{L}}+1$ equations of the
system (3)–(4). Since the indicial polynomial of $\mathcal{M}$ is a factor of
the indicial polynomial $p_{0}$, this system then implies that $T$ can be
extended to a unique power series solution of $\mathcal{M}$. As $\mathcal{M}$
is a right factor of $\mathcal{L}$, this power series is also a solution of
$\mathcal{L}$. Since it has the same initial conditions as $S$, they coincide.
∎

###  2.2. Homogeneous minimization

Since initial conditions are given for the power series $S$ solution of the
linear differential equation, it is possible to compute arbitrarily many
coefficients of $S$. The algorithm relies on the computation of upper bounds
on the degree of the coefficients of right factors of the linear differential
operator of a given order. Given such bounds and sufficiently many
coefficients of $S$, it is easy to set up a (structured) linear system whose
solutions are the possible coefficients of a right factor, or only 0 if no
such factor exists. When a non-zero solution is found, one takes its greatest
common right divisor with the original linear differential operator and checks
it using Lemma 2.1. This approach is described in Algorithm 1.

Algorithm 1 Minimal right factor

1:$\mathcal{L}=a_{r}(z)\partial_{z}^{r}+\dots+a_{0}(z)$ in
${\mathbb{Q}}[z]\langle\partial_{z}\rangle$; ini: $S_{0}$ a truncated power
series at precision $\geq\max\mathcal{Z}_{\mathcal{L}}$ specifying a unique
solution $S\in{\mathbb{Q}}[[z]]$ of $\mathcal{L}(S)=0$.

2:a right factor of $\mathcal{L}$ in
${\mathbb{Q}}[z]\langle\partial_{z}\rangle$ of minimal order that vanishes at
$S$

3:$\mathcal{M}:=\mathcal{L}$; $T:=S_{0}$; $m:=r$;
$p:=\max\mathcal{Z}_{\mathcal{L}}+r$

4:while $m>1$ do

5: $m:=m-1$

6: if $N:=$BoundDegreeCoeffs($\mathcal{L}$, $m$)$\neq$FAIL then

7: $p:=\max(p,m(N+1))$

8: while true do

9: $T:=$SeriesSolution($\mathcal{L},T,p+m$);

10:
$\mathcal{H}:=\textsc{ApproximantBasis}(T,T^{\prime},\dots,T^{(m)};N,\dots,N;p)$

11: if $\mathcal{H}=\emptyset$ then break//  _No right factor of order $m$_

12: else//  _$\mathcal{H}$ contains at least a candidate factor $h$_

13: $\mathcal{M}:=\textsc{GreatestCommonRightDivisor}(\mathcal{L},h)$;

14: if $\mathcal{M}(T)=O(z^{\max\mathcal{Z}_{\mathcal{L}}+1})$ then
$m:={\operatorname{ord}}\,\mathcal{M}$; break

15: else $p:=2p$

16:return $\mathcal{M}$

It relies on several other algorithms that we now review.

####  2.2.1. Sub-algorithms

##### SeriesSolution

Takes as input a linear differential operator, a truncated power series
solution of it, and a target precision $p$. It returns the power series
solution of the operator up to $O(z^{p})$, obtained either by truncating the
power series given as input, or by extending it using the linear recurrence
deduced from the differential equation.

##### ApproximantBasis

Takes as input ${k}$ power series $(S_{1},\dots,S_{k})$ that are the
truncations at precision $p$ of the successive derivatives of $S$; $k$
nonnegative integers $(s_{1},\dots,s_{k})$ and the precision $p$. It first
computes a basis $B(z)\in\mathbb{Q}[z]^{k\times k}$ of the
$\mathbb{Q}[z]$-module

$\mathcal{A}_{p}:=\\{(p_{1},\dots,p_{k})\mid
p_{1}S_{1}+\dots+p_{k}S_{k}=O(z^{p})\\}$ (5)

in _shifted Popov form_ [40, 47, 11, 29] with shift vector
$(-s_{1},\dots,-s_{k})$. This implies that any element $P$ of
$\mathcal{A}_{p}$ with degrees bounded by $(s_{1},\dots,s_{k})$ is a linear
combination of the rows of $B$ whose index $i$ satisfies $\deg B_{ii}\leq
s_{i}$. Those are the rows returned by ApproximantBasis. Efficient algorithms
to compute such bases are known [29].

##### GreatestCommonRightDivisor

Computes the monic greatest common right divisor (gcrd) of two linear
differential operators. This is classically achieved by a non-commutative
version of Euclid’s algorithm [37] and more efficient methods are known [27,
48].

##### BoundDegreeCoeffs

This is the heart of the algorithmic work, described in Section 2.3. It takes
as input an operator of order $r$ and a positive integer $m<r$. It returns
either FAIL when it has proved that no right factor of order $m$ with
polynomial coefficients exist; or an upper bound on the degree of each of the
coefficients such a factor would have.

######  Theorem 2.2.

Given a linear differential operator
$\mathcal{L}\in\mathbb{Q}[z]\langle\partial_{z}\rangle$ and a truncated power
series specifying a unique solution $S\in\mathbb{Q}[[z]]$ of
$\mathcal{L}(S)=0$, Algorithm 1 computes a non-zero right factor $\mathcal{M}$
of $\mathcal{L}$ of minimal order such that $\mathcal{M}(S)=0$.

###### Proof.

1\. (Correctness assuming termination.) Since $T$ is expanded at precision
$p+m$ in 9 and $p>\max\mathcal{Z}_{\mathcal{L}}$ from 3 and 7, it satisfies
$T^{(i)}(0)=S^{(i)}(0)$ for $i\in\mathcal{Z}_{\mathcal{L}}$. In 10, all series
$T,T^{\prime},\dots,T^{(m)}$ are known at precision $p$. It follows that if
the basis returned by ApproximantBasis is empty with the given bounds on the
degrees of the coefficients, there is no right-factor of $\mathcal{L}$ of
order $m$. Otherwise, taking $\mathcal{M}$ a gcrd of $\mathcal{L}$ and an
element $h$ of $\mathcal{H}$ gives a right factor of $\mathcal{L}$ to which
Lemma 2.1 applies, showing that $\mathcal{M}(S)=0$. The loop on $m$ makes the
algorithm terminate on a right factor of minimal order.

2\. (Termination.) The only possible source on non-termination in the
algorithm is the loop where $p$ is doubled every time $\mathcal{M}$ fails to
cancel $T$ to sufficient precision. Let $V_{p}$ be the $\mathbb{Q}$-vector
space generated by the approximants of the modules $\mathcal{A}_{p^{\prime}}$
from Eq. 5 for all $p^{\prime}\geq p$. Since the approximants have degrees
bounded by $(N,\dots,N)$, these are finite-dimensional vector spaces and
$V_{p+1}\subset V_{p}$. Thus there exists $p_{0}$ such that $V_{p_{0}}$ is the
intersection of all $V_{p}$ for $p\geq p_{0}$. Any approximant
$h=(h_{0},\dots,h_{m})$ in $\mathcal{H}$ in 10 for $p\geq p_{0}$ has the
property that $h_{0}S+\dots+h_{m}S^{(m)}=O(z^{k})$ for all $k\geq p$ and thus
annihilates $S$, and therefore so does its gcrd with $\mathcal{L}$, making the
algorithm terminate. ∎

####  2.2.2. Comparison with van Hoeij’s algorithm

Van Hoeij’s Algorithm “Construct $R$” [51, p. 552] follows a similar pattern.
Our termination proof is essentially his. The difference is that instead of
looking for an arbitrary right factor of $\mathcal{L}$, we need to make sure
that the factor returned by the algorithm cancels the power series $S$. This
is ensured by the test in 14.

####  2.2.3. Example

Consider the sequence

$u_{n}=\sum_{k=0}^{n}{\frac{n!(n+k)!}{k!^{4}(n-k)!^{3}}}.$

Zeilberger’s creative telescoping algorithm [54] shows that $u_{n}$ satisfies
a linear recurrence of order 4 with coefficients that are polynomials in $n$
of degree at most 10:

$\left(29412n^{4}+224352n^{3}+632931n^{2}+781692n+356309\right)\left(n+3\right)^{2}\left(n+4\right)^{4}u_{n+4}\\\
+\dots+4\left(29412n^{4}+342000n^{3}+1482459n^{2}+2838258n+2024696\right)\left(n+1\right)^{2}u_{n}=0.$

This recurrence translates into a linear differential operator of order $10$
annihilating the generating function $S(z)=\sum_{n\geq 0}{u_{n}z^{n}}$, with
coefficients of degree at most 8:

$\mathcal{L}=29412z^{8}\partial_{z}^{10}-684z^{7}\left(688z-1489\right)\partial_{z}^{9}-21z^{6}\left(156864z^{2}+742368z-588707\right)\partial_{z}^{8}+\dotsb\\\
+\left(99370416z^{3}-1926228512z^{2}-19342508z+8500\right)\partial_{z}+4\left(2024696z^{2}-3141504z-32725\right).$
(6)

The only integer roots of the indicial polynomial of $\mathcal{L}$ at 0 are in
$\mathcal{Z}_{\mathcal{L}}=\\{0,1\\}$ so that the initial conditions
specifying $S$ uniquely are $S(0)=u_{0}=1,S^{\prime}(0)=u_{1}=3$. The
differential operator $\mathcal{L}$ is not minimal for $S$. There are two
stages in the execution of the algorithm: first, a right factor is sought;
next, its minimality is proved.

The algorithm first tries to find a right factor of $\mathcal{L}$ order 9 and
discovers while computing bounds on the degrees of its coefficients that no
such factor exists. The same observation is made for orders 8 and 7. For order
$6$, a bound 30 for the degrees of the coefficients is found. With this bound,
a candidate linear differential operator of order 6 is found:

$\mathcal{M}=z^{4}\left(1882368z^{4}-2206584z^{3}+1703460z^{2}+67815z+272\right)\partial_{z}^{6}+\dotsb+\\\
+2\left(3764736z^{6}-41001696z^{5}+157022376z^{4}-184937064z^{3}-6917519z^{2}-3408891z-41888\right).$

The computation of the greatest common right factor stops at its first step,
discovering that this differential operator is a right factor of
$\mathcal{L}$.

At order 5, a bound $N=15$ on the degree of the coefficients of possible
factors is computed. With $p+5$ coefficients of $S$, where $p=80=5\times 16$,
the computation of ApproximantBasis shows that there is no non-zero operator
$h$ of order 5 with coefficients of degree at most $N$ such that
$h(S(z))=O(z^{p})$ and therefore no right factor of $\mathcal{L}$ of order 5
annihilating $S$. Next, the bound on the degrees of the coefficients of a
right factor of order $4$ is smaller than $15$, so that if a right factor of
that order existed, it would have been obtained for order 5. Finally, the
computations of bounds for orders $3,2,1$ show that no factor of these orders
exist. This concludes the proof of minimality of the operator $\mathcal{M}$
for $S$.

###  2.3. Degree bounds

The computation of degree bounds for a factor of a given order is a key step
in van Hoeij’s factorization algorithm [50, §9]. We recall the ingredients
here. Compared to our earlier work [17] where we have obtained universal
bounds, the bounds computed here are tailored to the equation under study,
rather than depending only on its order, degree and height. This allows for
smaller bounds and more efficient computations.

####  2.3.1. Singularities of the factors

Dividing $\mathcal{L}$ by its leading coefficient $a_{r}$ gives a monic
operator with rational function coefficients. In this form, the singularities
of $\mathcal{L}$ are the poles of its coefficients. A singularity $\alpha$ of
$\mathcal{L}$ is called _regular_ if the indicial polynomial of $\mathcal{L}$
at $\alpha$ has degree equal to the order $r$ of $\mathcal{L}$, and it is
called _irregular_ otherwise. The right factors will be searched in the same
monic form. Recall that the valuation ${\operatorname{val}}_{\alpha}(r)$ of a
rational function $r$ at $\alpha$ is the exponent of the leading term of the
Laurent expansion of $r$ at $\alpha$ (and
${\operatorname{val}}_{\alpha}(0)=\infty)$. At a regular singularity the
valuation of each coefficient $a_{i}$ of $\mathcal{L}$ in monic form is at
least $i-r$. Bounds on the degrees of the coefficients of factors are obtained
by bounding the valuations of their coefficients in monic form at each
singularity and at infinity, and by bounding the number of _apparent_
singularities. Apparent singularities are poles of the coefficients where the
operator has a basis of $r$ formal power series solutions; they are regular.
All these notions are classical and can be found for instance in Ince’s book
[28].

####  2.3.2. Newton polygons and valuations of the coefficients of the factors

The Newton polygon of the operator $\mathcal{L}$ from Eq. 1 at 0 is the convex
hull of the union of the quadrants
$(i,{\operatorname{val}}_{0}(a_{i})-i)+(\mathbb{R}_{\leq
0}\times\mathbb{R}_{\geq 0})$. The knowledge of the Newton polygon of
$\mathcal{L}$ at 0 gives lower bounds on the valuations of its coefficients.
The main property of relevance here is that the Newton polygon of a product of
operators is the (Minkowski) sum of their Newton polygons ([33, Lemme 1.4.1]).
For instance, when 0 is an ordinary point or a regular singularity of
$\mathcal{L}$, the only slope of the lower part of its Newton polygon is 0 and
this is therefore a property of the Newton polygons of the monic factors of
$\mathcal{L}$, which reflects the fact that they are regular at 0 in that
case.

More generally, let $(x_{0},y_{0})=(0,y_{0}),\dots,(x_{k},y_{k})=(r,y_{k})$ be
the points on the lower part of the Newton polygon of $\mathcal{L}$ and let
$((n_{1},d_{1}),\dots,(n_{k},d_{k}))$ with
$(n_{i},d_{i})=(x_{i}-x_{i-1},y_{i}-y_{i-1})$ be the tuple of segments of the
Newton polygon of $\mathcal{L}$ sorted by increasing slope. Then the lowest
possible Newton polygon for a monic factor of order $m$ is obtained from the
solution of the “0-1 knapsack problem”

$\min\sum_{i=1}^{k}{c_{i}n_{i}}\quad\text{subject to
}\sum_{i=1}^{k}{c_{i}d_{i}}=m\text{ and }c_{i}\in\\{0,1\\},\quad i=1,\dots,k,$

where $c_{i}$ is either 1 or 0 depending on whether or not the slope
$(n_{i},d_{i})$ is used. This solution allows one to obtain lower bounds on
the valuations at 0 of the coefficients of monic factors of $\mathcal{L}$ of
order $m$. The 0-1 knapsack problem is NP-hard but lower bounds can be found
efficiently if needed [52, Ch. 8]. In practice, this has never been a costly
step in our computations and an optimal value may lead to a better degree
bound which saves computation time in other steps of the minimization
algorithm.

The same process can be performed at every irregular singularity $\alpha$ of
$\mathcal{L}$ by considering the Newton polygon formed from
${\operatorname{val}}_{\alpha}$ instead of ${\operatorname{val}}_{0}$. Thus
lower bounds on the valuations of the coefficients of a factor are found at
each singularity, from the Newton polygon of $\mathcal{L}$ and the order of
the factor. Applying the same process at $\infty$ (for instance by changing
$z$ into $1/z$ and working at 0) gives bounds on the valuation at infinity of
these coefficients.

####  2.3.3. Fuchs relation and apparent singularities of the factors

Let $\mathcal{M}$ denote a monic right-factor of order $m$ of the operator
$\mathcal{L}$ to be minimized. The study of the Newton polygons of
$\mathcal{L}$ provides lower bounds on the valuations of the coefficients of
$\mathcal{M}$ at the singularities of $\mathcal{L}$. We now show how to obtain
an upper bound on the number of apparent singularities of $\mathcal{M}$;
together with the lower bounds on valuations, this will provide upper bounds
on the degrees of (the polynomial version of) $\mathcal{M}$.

We first recall the principle of the method in the case where $\mathcal{M}$ is
Fuchsian, that is, if all its singularities (including $\infty$) are regular.
Fuchs’ relation [49, p. 138] states that

$\sum_{\rho\in\operatorname{Sing}(\mathcal{M})}S_{\rho}(\mathcal{M})=-m(m-1),$
(7)

where $\operatorname{Sing}(\mathcal{M})$ is the set of singularities of
$\mathcal{M}$, including the apparent ones and infinity, and where

$S_{\rho}(\mathcal{M}):=\sum_{j=1}^{m}{e_{j}(\rho)}-\frac{m(m-1)}{2},$ (8)

the numbers $e_{j}(\rho)$ being the local exponents of $\mathcal{M}$ at the
point $\rho$ (they are the roots of the indicial polynomial at $\rho$). At an
apparent singularity $\rho$, the quantity $S_{\rho}(\mathcal{M})$ is a
positive integer, so that the number of apparent singularities is upper
bounded by

$\\#\\!\operatorname{App}(\mathcal{M})\leq-m(m-1)-\sum_{\rho\in\sigma(\mathcal{M})}S_{\rho}(\mathcal{M}),$
(9)

with $\sigma(\mathcal{M})$ the subset of $\operatorname{Sing}(\mathcal{M})$
formed by the singularities of $\mathcal{M}$ that are not apparent.

Since $\mathcal{M}$ is a right factor of $\mathcal{L}$, the set
$\sigma(\mathcal{M})$ is a subset of $\operatorname{\sigma(\mathcal{L})}$. The
set $\operatorname{\sigma(\mathcal{L})}$ is known, since it corresponds to the
roots of the leading coefficient $a_{r}(z)$ that are not apparent
singularities of $\mathcal{L}$, plus possibly $\infty$. Let
$\mu_{1},\dots,\mu_{s}$ be the irreducible factors of $a_{r}(z)$ corresponding
to non-apparent singularities of $\mathcal{L}$ and by convention let
$\mu_{0}=z$. At a finite $\rho\in\sigma(\mathcal{L})$, given by its minimal
polynomial $\mu_{i}$, the indicial polynomial
$\operatorname{ind}_{\rho}^{\mathcal{L}}(\theta)\in{\mathbb{Q}}(\rho)[\theta]$
is easily computed. Then the unknown indicial polynomial
$\operatorname{ind}_{\rho}^{\mathcal{M}}(\theta)\in{\mathbb{Q}}(\rho)[\theta]$
has to be a factor of $\operatorname{ind}_{\rho}^{\mathcal{L}}(\theta)$ of
degree exactly $m$. Let $I_{i,j}(\theta)$, $j=1,\dots,k_{i}\leq r$ be the
irreducible factors of $\operatorname{ind}_{\rho}^{\mathcal{L}}(\theta)$ in
${\mathbb{Q}}(\rho)[\theta]$, repeated with their multiplicity (and similarly
$I_{0,j}(\theta)$ denote the factors of the indicial polynomial of
$\mathcal{L}$ at infinity). The sum of the roots of $I_{i,j}$ lies in
${\mathbb{Q}}(\rho)$ and its sum over all roots of $\mu_{i}$ is a rational
number $e_{i,j}$. A bound on the number of apparent singularities is therefore
obtained by solving the following 0-1 knapsack problem

$\displaystyle\max A\quad\text{subject to }\
A=-m(m-1)-\sum_{i=0}^{s}\deg\mu_{i}\left(\sum_{j=1}^{k_{i}}{c_{i,j}e_{i,j}}-\frac{m(m-1)}{2}\right)\in\mathbb{N}$
$\displaystyle\text{and for all
$i\in\\{0,\dots,s\\}$,}\quad\sum_{j=1}^{k_{i}}{c_{i,j}\deg I_{i,j}}=m,\quad
c_{i,j}\in\\{0,1\\}.$

The constraints express the fact that there should be $m$ exponents at each
root of $a_{r}$, be them singular or ordinary for $\mathcal{M}$.

Note that if there is no choice of $c_{ij}$ for which $A\in\mathbb{N}$, then
there is no right factor of order $m$. Also, as for the previous optimization
problem, this is potentially a computationally expensive step and simple upper
bounds can be obtained efficiently, but this step did not seem to be an
obstacle in our experiments.

This process can be used whenever the right factor $\mathcal{M}$ to be found
is known to be Fuchsian, thus in particular when $\mathcal{L}$ itself is
Fuchsian.

####  2.3.4. Generalized Fuchs relation

To an irregular singular point $\rho$ of $\mathcal{L}$ is associated a set of
_exponential parts_. If $\rho$ is finite, these are polynomials $w(z)$ in some
rational power $1/r$ of $z$ ($r\in\mathbb{N}_{>0}$) such that $\mathcal{L}$
admits a formal solution of the form

$\exp\\!\left(\int{\frac{w(1/(z-\rho))}{z-\rho}\,dz}\right)S(z),\qquad
S\in{\mathbb{Q}}\\!\left[\left[(z-\rho)^{1/r}\right]\right][\log(z-\rho)],\quad\operatorname{val}_{z=\rho}S=0.$

The case when $\rho=\infty$ is obtained by changing $z$ into $1/z$ in the
equation. The computation of the list of exponential parts at a point $\rho$
is well understood [33]. When $\rho$ is a regular singular point, then the
exponential parts are constants that coincide with the roots of the indicial
polynomial. In general, the _generalized local exponents_ are the constant
coefficients of the exponential parts. If each of them is counted with
multiplicity $r$, their number is exactly the order of the operator.

When $\mathcal{M}$ is a right factor of $\mathcal{L}$, its exponential parts
at $\rho$ form a subset of those of $\mathcal{L}$. If the order of
$\mathcal{M}$ is $m$, the Fuchs relation (7) generalizes as

$\sum_{\rho\in\operatorname{Sing}(\mathcal{M})}\\!\left(S_{\rho}(\mathcal{M})-\frac{1}{2}I_{\rho}(\mathcal{M})\right)=-m(m-1),$
(10)

where $S_{\rho}$ is as in Eq. 8, with the generalized local exponents taking
the place of the local exponents and

$I_{\rho}(\mathcal{M}):=2\sum_{1\leq i<j\leq m}\deg(w_{i}-w_{j}),$

where the $w_{i}$ are the exponential parts at $\rho$, see [12]. As the
$w_{i}$ are polynomials in a fractional power of $z$, their degree here is a
rational number. (The quantity $I_{\rho}(\mathcal{M})$ is related to
Malgrange’s irregularity of $\mathcal{M}$ at $\rho$; see [17, §2.2] for
details which are not essential here.)

Thus, the analogue of Eq. 9 is

$\\#\\!\operatorname{App}(\mathcal{M})\leq-m(m-1)-\sum_{\rho\in\sigma(\mathcal{M})}\left(S_{\rho}(\mathcal{M})-\frac{1}{2}I_{\rho}(\mathcal{M})\right).$
(11)

The corresponding optimization problem is slightly more involved. As in the
Fuchsian situation, $\sigma(\mathcal{M})\subset\sigma(\mathcal{L})$ and we
denote by $\mu_{1},\dots,\mu_{s}$ the irreducible factors of the leading
coefficient $a_{r}(z)$ that correspond to non-apparent singularities of
$\mathcal{L}$ and $\mu_{0}=z$. At a finite $\rho\in\sigma(\mathcal{L})$, given
by its minimal polynomial $\mu_{i}$, the exponential parts are given as

$w_{i1}((z-\rho)^{1/r_{i1}}),\dots,w_{ik_{i}}((z-\rho)^{1/r_{ik_{i}}})$

with minimal $r_{ij}$. Each contributes $r_{ij}$ times its constant
coefficient to the set of generalized local exponents of $\mathcal{L}$ at
$\rho$, so that
$\sum{r_{ij}}=\operatorname{{\operatorname{ord}}}(\mathcal{L}).$ The
exponential parts of $\mathcal{M}$ at $\rho$ form a subset of those of
$\mathcal{L}$. This property, combined with the generalized Fuchs relation
(10), leads to the following linear optimization problem

$\displaystyle\max A\quad{\text{subject to}}$
$\displaystyle\begin{split}A=-m(m-1)-\sum_{i=0}^{s}\sum_{\rho\mid\mu_{i}(\rho)=0}\left(\sum_{j=1}^{k_{i}}{{c_{i,j}(r_{i,j}w_{i,j}(0)-\frac{r_{ij}(r_{ij}-1)}{2}\deg
w_{ij})-\frac{m(m-1)}{2}}}\right)\\\
+\sum_{i=0}^{s}\deg\mu_{i}\sum_{j=1}^{k_{i}}\sum_{1\leq k\neq j\leq
k_{i}}d_{i,\\{j,k\\}}\deg(w_{ij}-w_{ik})\in\mathbb{N}\end{split}$
$\displaystyle\text{and for all
$i\in\\{0,\dots,s\\}$,}\quad\sum_{j=1}^{k_{i}}{c_{i,j}r_{i,j}}=m,\quad
c_{i,j}\in\\{0,1\\},$ $\displaystyle\text{and for all
$(i,j)$,}\quad\sum_{k\neq j}d_{i,\\{j,k\\}}=c_{i,j}(m-1),\quad
d_{i,\\{j,k\\}}\in\\{0,1\\}.$

The last set of constraints consists in adding one variable for each pair of
$(w_{ij},w_{ik})$ and forcing the sum of these variables for fixed $i$ to be
the number of pairs, namely $m-1$, an idea taken from [20].

####  2.3.5. Example

Consider the equation

$zy^{\prime\prime}+(1-6z)y^{\prime}+(z-3)y=0,$

with initial condition $y(0)=1$, which specifies a unique power series
solution $S(z)=1+3z+13z^{2}/2+\dotsb$. It has two singular points, at 0 and
$\infty$. The point 0 is regular, with exponents $0,0$. The point $\infty$ is
irregular, with exponential parts $w_{\pm}=\alpha_{\pm}z+1/2$, where
$\alpha_{\pm}=-3\pm 2\sqrt{2}$, corresponding to formal solutions
$\exp(-\alpha_{\pm}z)/\sqrt{z}$ at infinity and both generalized exponents are
equal to $1/2$. In the notation above, we have

$\displaystyle s=1,\quad\mu_{0}=\mu_{1}=z,\quad k_{0}=k_{1}=2,\quad
r_{0,1}=r_{0,2}=r_{1,1}=r_{1,2}=1,$ $\displaystyle
w_{0,1}=w_{+},w_{0,2}=w_{-},w_{1,1}=w_{1,2}=0.$

Looking for a right factor of order $m$ leads to maximizing $A\in\mathbb{N}$
such that

$A=-m(m-1)-\left(c_{0,1}/2+c_{0,2}/2-m(m-1)/2\right)+d_{0,\\{1,2\\}},$

with the constraints

$\displaystyle c_{0,1}+c_{0,2}=m,\quad
d_{0,\\{1,2\\}}=c_{0,1}(m-1)=c_{0,2}(m-1),$ $\displaystyle
c_{0,1},c_{0,2},d_{0,\\{1,2\\}}\text{ in }\\{0,1\\}.$

For the order $m=1$ of a right factor, the last constraints force
$d_{0,\\{1,2\\}}=0$, $c_{0,1}+c_{0,2}=1$, which makes $A<0$, showing that
there is no solution and thus no factorization with a right factor of order 1;
the equation is minimal.

####  2.3.6. Example

We show in more detail the computation for the order 10 differential equation
(6) of Section 2.2.3.

There are two singularities: 0 and infinity. The point 0 is regular, with
exponents

$0,0,0,0,1,1,\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},$

with $\alpha_{i}$ the four roots of the irreducible polynomial

$P_{\alpha}=29412x^{4}-246240x^{3}+764259x^{2}-1042332x+527381.$

The point $\infty$ is _irregular_. Its exponential parts are

$1,1,\beta_{i}\ (i=1,\dots,4),\gamma_{i}x+3/2\ (i=1,\dots,4),$

with $\beta_{i}$ and $\gamma_{i}$ roots of the irreducible polynomials

$P_{\beta}=29412x^{4}-342000x^{3}+1482459x^{2}-2838258x+2024696,\quad
P_{\gamma}=x^{4}+16x^{3}-112x^{2}+284x+4.$

Thus, for this equation,

$S_{0}(\mathcal{L})=2+\sum_{i}{\alpha_{i}}-45=-\frac{1489}{43},\quad
S_{\infty}(\mathcal{L})=2+\sum_{i}{\beta_{i}}+4\frac{3}{2}-45=-\frac{1901}{43},\quad
I_{\infty}(\mathcal{L})=60$

and the generalized Fuchs equation reduces to

$-\frac{1489}{43}-\frac{1901}{43}-30=-90.$

For lower orders $m$, the optimization problem to be solved is

$\displaystyle\max A\quad{\text{subject to}}$
$\displaystyle\begin{split}A=-m(m-1)-\left(c_{0,1}+c_{0,2}+c_{0,\\{3,4,5,6\\}}\frac{500}{43}+c_{0,\\{7,8,9,10\\}}4\frac{3}{2}-\frac{m(m-1)}{2}\right)\\\
-\left(c_{1,5}+c_{1,6}+c_{1,\\{7,8,9,10\\}}\frac{360}{43}-\frac{m(m-1)}{2}\right)\\\
+4d_{0,\\{1,\\{7,8,9,10\\}\\}}+4d_{0,\\{2,\\{7,8,9,10\\}\\}}+16d_{0,\\{\\{3,4,5,6\\},\\{7,8,9,10\\}\\}}+6d_{0,\\{\\{7,8,9,10\\},\\{7,8,9,10\\}\\}}\in\mathbb{N}\end{split}$
with constraints $\displaystyle
c_{0,1}+c_{0,2}+4c_{0,\\{3,4,5,6\\}}+4c_{0,\\{7,8,9,10\\}}=c_{1,1}+c_{1,2}+c_{1,3}+c_{1,4}+c_{1,5}+c_{1,6}+4c_{1,\\{7,8,9,10\\}}=m$
$\displaystyle
d_{0,\\{1,2\\}}+4d_{0,\\{1,\\{3,4,5,6\\}\\}}+4d_{0,\\{1,\\{7,8,9,10\\}\\}}=c_{0,1}(m-1),$
$\displaystyle
d_{0,\\{1,2\\}}+4d_{0,\\{2,\\{3,4,5,6\\}\\}}+4d_{0,\\{2,\\{7,8,9,10\\}\\}}=c_{0,2}(m-1),$
$\displaystyle
d_{0,\\{1,\\{3,4,5,6\\}\\}}+d_{0,\\{2,\\{3,4,5,6\\}\\}}+3d_{0,\\{3,4,5,6\\},\\{3,4,5,6\\}}+4d_{0,\\{\\{3,4,5,6\\},\\{7,8,9,10\\}\\}}=(m-1)c_{0,\\{3,4,5,6\\}},$
$\displaystyle
d_{0,\\{1,\\{7,8,9,10\\}\\}}+d_{0,\\{2,\\{7,8,9,10\\}\\}}+4d_{0,\\{\\{3,4,5,6\\},\\{7,8,9,10\\}\\}}+3d_{0,\\{7,8,9,10\\},\\{7,8,9,10\\}}=(m-1)c_{0,\\{7,8,9,10\\}},$
$\displaystyle\text{and for all $(i,j,k)$,}\qquad c_{i,j}\in\\{0,1\\},\quad
d_{i,\\{j,k\\}}\in\\{0,1\\}.$

Integrality of $A$ forces $c_{0,\\{3,4,5,6\\}}=c_{1,\\{7,8,9,10\\}}$. If they
are both equal to 1, the first two lines of the constraint on $A$ give a
quantity that is at most $-20$. Making $A\geq 0$ then requires
$d_{0,\\{3,4,5,6\\},\\{7,8,9,10\\}}=1$. The last constraint then makes
$c_{0,\\{7,8,9,10\\}}=1$, which turns the constraint on $A$ into

$-20-6+16+4d_{0,\\{1,\\{7,8,9,10\\}\\}}+4d_{0,\\{2,\\{7,8,9,10\\}\\}}+6d_{0,\\{\\{7,8,9,10\\},\\{7,8,9,10\\}\\}}\geq
0.$

Therefore $d_{0,\\{\\{7,8,9,10\\},\\{7,8,9,10\\}\\}}=1$ and at least one of
$d_{0,\\{1,\\{7,8,9,10\\}\\}}$ and $d_{0,\\{2,\\{7,8,9,10\\}\\}}$ is 1 too.
The last constraint then shows that $m=9$ or $m=10$ depending on whether one
or two of them are 1. We know that $m=10$ is possible: it is the original
equation. If $m=9$, then the first constraint gives $c_{0,1}+c_{0,2}=1$.
Injecting into the constraint for $A$ makes $A<0$, a contradiction.

We have thus proved that for a strict factor of $\mathcal{A}$,
$c_{0,\\{3,4,5,6\\}}=c_{1,\\{7,8,9,10\\}}=0$. This makes all variables 0 in
the left-hand side of the penultimate constraint. The last constraint becomes

$d_{0,\\{1,\\{7,8,9,10\\}\\}}+d_{0,\\{2,\\{7,8,9,10\\}\\}}+3d_{0,\\{7,8,9,10\\},\\{7,8,9,10\\}}=(m-1)c_{0,\\{7,8,9,10\\}}.$

If $c_{0,\\{7,8,9,10\\}}$ was equal to 0, then the constraint on $A$ would
give $c_{0,1}=c_{0,2}=0$ too, which would give $m=0$ in the second one, a
contradiction. Therefore $c_{0,\\{7,8,9,10\\}}=1$. The remaining constraints
are

$\displaystyle\begin{split}A=-\left(c_{0,1}+c_{0,2}+6+c_{1,5}+c_{1,6}\right)\qquad\qquad\qquad\qquad\qquad\qquad\qquad\\\
+4d_{0,\\{1,\\{7,8,9,10\\}\\}}+4d_{0,\\{2,\\{7,8,9,10\\}\\}}+6d_{0,\\{\\{7,8,9,10\\},\\{7,8,9,10\\}\\}}\geq
0,\end{split}$ $\displaystyle
c_{0,1}+c_{0,2}+4=c_{1,1}+c_{1,2}+c_{1,3}+c_{1,4}+c_{1,5}+c_{1,6}=m,$
$\displaystyle d_{0,\\{1,2\\}}+4d_{0,\\{1,\\{7,8,9,10\\}\\}}=c_{0,1}(m-1),$
$\displaystyle d_{0,\\{1,2\\}}+4d_{0,\\{2,\\{7,8,9,10\\}\\}}=c_{0,2}(m-1),$
$\displaystyle
d_{0,\\{1,\\{7,8,9,10\\}\\}}+d_{0,\\{2,\\{7,8,9,10\\}\\}}+3d_{0,\\{7,8,9,10\\},\\{7,8,9,10\\}}=m-1.$

The second one then implies that _the order of a strict right factor of
$\mathcal{L}$ can only be one of $\\{4,5,6\\}$._

With $m=6$, there is only one solution (meaning that no optimization is
needed), with all the remaining variables equal to 1 and the bound $A$ on the
number of apparent singularities equal to 4. There are therefore at most 5
regular singularities (these four and 0, which is a regular singularity of
$\mathcal{L}$). Such a factor can be written

$\mathcal{M}=\partial_{z}^{6}+\frac{a_{5}}{A}\partial_{z}^{5}+\dots+\frac{a_{0}}{A^{6}},$

with $A$ of degree at most 5. The Newton polygon of $\mathcal{L}$ at infinity
has for vertices $(0,0),(6,0),(10,4)$. The largest possibility for
$\mathcal{M}$ is therefore $(0,0),(2,0),(6,4)$, leading to the following
bounds on the degree of the numerators $a_{i}$: $(\deg A-1,\deg A^{2}-2,\deg
A^{3}-4,\deg A^{4}-4,\deg A^{5}-4,\deg A^{6}-4)$. Reducing to the same
denominator gives the bounds $(30,29,28,26,26,26,26)$ on the degrees of the
coefficients of $(\partial^{6},\dots,1)$. This is the bound used in Example
2.2.3, leading to the discovery of the factor $\mathcal{M}$.

With $m=5$, there are several solutions, which are as in the case when $m=6$,
but with one of $c_{0,1}$ and $c_{0,2}$ equal to 0 and consequently
$d_{0,\\{1,2\\}}=0$ and one of $d_{0,\\{1,\\{7,8,9,10\\}\\}}$,
$d_{0,\\{2,\\{7,8,9,10\\}\\}}$ equals 0, leading to a bound $A\leq
6+4-(6+1+1)=2$ on the number of apparent singularities. By the same reasoning
as above, this leads to the bounds $(15,14,13,13,13,13)$ on the degrees of the
coefficients of a factor of order 5. Using about 90 coefficients of the series
shows that such a factor does not exist.

Finally, with $m=4$ the only of the remaining $d$ variables that is not 0 is
$d_{0,\\{7,8,9,10\\},\\{7,8,9,10\\}}$ and the bound on $A$ becomes $6-(6)=0$.
Again, a computation with degree bounds $(5,4,3,3,3)$ proves that no factor of
degree 4 exists.

###  2.4. Inhomogeneous minimization

Again, we consider an equation like Eq. 1 and initial conditions for a unique
formal power series $S$ solution of it. Using the method of the previous
section, we can assume that it has minimal order. The problem of inhomogeneous
minimisation is to find an equation

$\mathcal{M}(y(z))=B(z),\quad\text{with}\quad\mathcal{M}=b_{s}(z)\partial_{z}^{s}+\dots+b_{0}(z),$

with $s<r$ and rational function coefficients $b_{0},\dots,b_{s},B$
($b_{s}\neq 0$), having $S$ as a solution. When such an equation exists with
$B\neq 0$, applying $B\partial_{z}-B^{\prime}$ on both sides of the equation
yields a homogeneous linear differential equation of order $s+1$ satisfied by
$S$, so that minimality of $\mathcal{L}$ implies $s=r-1$. Without loss of
generality (up to replacing $\mathcal{M}$ by $\frac{1}{B}\mathcal{M}$) one can
assume $B(z)=1$ and then differentiation implies

$\partial\mathcal{M}=R(z)\mathcal{L}$

for some non-zero rational function $R$, which is therefore an _integrating
factor_ of $\mathcal{L}$. This implies that $R$ is a rational function
solution of the adjoint equation [39, Chap. III.§10]

$\mathcal{L}^{*}(R)=0.$ (12)

Finding rational solutions of linear differential equations is a classical
problem, whose solution can be found by an algorithm due to Abramov [1, 2],
with roots in Liouville’s work [31]. This algorithm returns a basis of
rational solutions of Eq. 12. This is a decision algorithm: if no non-zero
rational solution is found, this proves that there is no inhomogeneous linear
differential equation of order smaller than $r$ satisfied by the power series
$S$. Otherwise, minimality implies that the basis consists of one solution
$R(z)$. The operator $\mathcal{M}$ (known as the _bilinear concomitant_ [39])
can be reconstructed coefficient by coefficient (this is equivalent to [5, p.
703]). Then by design, $\mathcal{M}(S)$ is a constant $c$, which can be
computed from the first coefficients of the power series $S$, completing the
computation of the minimal inhomogeneous equation $\mathcal{M}(y)=c$ satisfied
by $S$. This computation is summarized in Algorithm 2.

Algorithm 2 Minimal inhomogeneous linear differential equation

1: $\mathcal{L}=a_{r}(z)\partial_{z}^{r}+\dots+a_{0}(z)$, a linear operator of
minimal order that vanishes at $S$; ini: $S_{0}$ a truncated power series at
precision $\geq\max\mathcal{Z}_{\mathcal{L}}$

2: $\mathcal{M}=b_{r-1}(z)\partial_{z}^{r-1}+\dots+b_{0}(z)$ and
$B(z)\in\mathbb{Q}(z)$, such that $\mathcal{M}(S)=B$ or Fail if no such pair
exists.

3:$\mathcal{L}^{\star}:=\operatorname{adjoint}(\mathcal{L})$;

4:$\mathcal{S}=\textsc{BasisRationalSolutions}(\mathcal{L}^{\star})$

5:if $\mathcal{S}=\emptyset$ then return FAIL

6:Let $R$ be the unique element of $\mathcal{S}$

7:$b_{r-1}:=Ra_{r}$

8:for $j=r-2,\dots,0$ do $b_{j}:=Ra_{j+1}-b_{j+1}^{\prime}$

9:Compute $S$ up to precision $r-\min_{j}\operatorname{val}_{0}(b_{j})$

10:Let $B$ be the constant term of $\mathcal{M}(S)$

11:return $\mathcal{M},B$

##  3\. Efficient computation of the set of algebraic values taken by
$E$-functions at algebraic points

###  3.1. The Adamczewski-Rivoal algorithm

In this section, we consider mainly $E$-functions with Taylor coefficients in
$\mathbb{Q}$. An _$E$ -function over $\mathbb{Q}$_ is a power series

$f(z):=\sum_{n=0}^{\infty}\frac{a_{n}}{n!}z^{n}\text{ in }{\mathbb{Q}}[[z]]$

such that there exists $C>0$ with the following properties:

  1. $(i)$

$f$ satisfies a homogeneous linear differential equation with coefficients in
${\mathbb{Q}}(z)$;

  2. $(ii)$

for any $n\geq 0$, $|a_{n}|\leq C^{n+1}$;

  3. $(iii)$

for any $n\geq 0$, there exists $d_{n}\in\mathbb{N}\setminus\\{0\\}$ such that
$d_{n}\leq C^{n+1}$ and $d_{n}a_{m}\in{\mathbb{Z}}$ for all $0\leq m\leq n$.

We shall sometimes simply write “ODE” for “linear differential equation with
coefficients in ${\mathbb{Q}}(z)$ or in ${\overline{\mathbb{Q}}}(z)$”; unless
otherwise stated, an ODE will be assumed to be homogeneous. In the rest of
this section, $E$-functions over $\mathbb{Q}$ are simply called $E$-functions.
This is justified by the fact that most of the discussion applies to more
general settings, in particular to $E$-functions with Taylor coefficients in
${\overline{\mathbb{Q}}}$ and to $E$-functions in Siegel’s more general sense,
as discussed in Sections 3.4 and 3.5.

$E$-functions are entire functions (by $(ii)$). Polynomials in
${\mathbb{Q}}[z]$ are trivial examples of $E$-functions; all non-polynomial
$E$-functions are transcendental over ${\mathbb{Q}}(z)$. The class of
$E$-functions includes the exponential function $\exp(z)$, Bessel’s function
of the first kind

$J_{0}(z):=\sum_{m=0}^{\infty}\frac{\left(-1\right)^{m}}{m!^{2}}\left(\frac{z}{2}\right)^{2m}=\,_{0}F_{1}[\,\cdot\,;1;-z^{2}/4],$

and more generally the hypergeometric $E$-functions, i.e. series of the form

${}_{p}F_{q}[a_{1},\ldots,a_{p};b_{1},\ldots,b_{q};\lambda
z^{q-p+1}]:=\sum_{n=0}^{\infty}\frac{(a_{1})_{n}\cdots(a_{p})_{n}}{(1)_{n}(b_{1})_{n}\cdots(b_{q})_{n}}\lambda^{n}z^{(q-p+1)n}$

with rational parameters $a_{i},b_{j}$, $q\geq p\geq 0$,
$\lambda\in{\overline{\mathbb{Q}}}^{*}$ and where $(\alpha)_{0}:=1$,
$(\alpha)_{n}:=\alpha(\alpha+1)\cdots(\alpha+n-1)$ for $n\geq 1$.
$E$-functions form a sub-ring of the ring of formal power series in
${\mathbb{Q}}[[z]]$, stable by $d/dz$ and $\int_{0}^{z}$; these properties can
be used to construct many examples of $E$-functions starting from
hypergeometric series. Shidlovskii has proved in [43, p. 184] that any
$E$-function solution of an ODE of order $1$ is of the form $p(z)e^{\lambda
z}$ for some $p(z)\in{\overline{\mathbb{Q}}}[z]$ and
$\lambda\in{\overline{\mathbb{Q}}}$. Gorelov has proved in [26] that any
$E$-function solution of an ODE of order $2$ is a
${\overline{\mathbb{Q}}}(z)$-linear combination of hypergeometric
$E$-functions with $p=q=1$ (he had obtained earlier in [25] a similar but more
precise result for $E$-functions solution of an inhomogeneous ODE of order
$1$). However, Fresán and Jossen have recently showed in [24] that not all
$E$-functions are ${\overline{\mathbb{Q}}}(z)$-linear combinations of
hypergeometric $E$-functions, nor even more generally polynomials in
hypergeometric $E$-functions with algebraic coefficients.

As of today, no algorithm is known neither for deciding whether a linear
differential equation $\mathcal{L}(y(z))=0$ admits solutions that are
$E$-functions, nor for deciding whether a solution of $y(z)$ of $\mathcal{L}$,
uniquely determined by sufficiently many initial conditions, is an
$E$-function. It is actually not clear whether these questions are decidable
or not. Consequently, the algorithm described below relies on the following
assumption:

  1. (A)

An oracle guarantees that the input $f$ is an $E$-function.

In practice, an $E$-function is given by an explicit expression for its Taylor
coefficients as a multiple hypergeometric sum and $\mathcal{L}$ can then be
computed for instance by Zeilberger’s creative telescoping algorithm [54].

Siegel initiated a program to determine when an $E$-function takes a
transcendental value at an algebraic point [44]. This culminated with the
celebrated Siegel-Shidlovskii theorem: given a vector $Y$ of $E$-functions
$f_{1},\ldots,f_{n}$ solution of a differential system $Y^{\prime}=AY$ with a
matrix $A$ with elements in ${\overline{\mathbb{Q}}}(z)$, the transcendence
degree over ${\overline{\mathbb{Q}}}(z)$ of the field generated by
$f_{1}(z),\ldots,f_{n}(z)$ over ${\overline{\mathbb{Q}}}(z)$ is equal to the
transcendence degree over ${\overline{\mathbb{Q}}}$ of the field generated by
$f_{1}(\alpha),\ldots,f_{n}(\alpha)$ over ${\overline{\mathbb{Q}}}$ for every
non-zero algebraic number $\alpha$ which is not a singularity of $A$ (i.e., a
pole of some element of $A$). In 2006, Beukers [14, Thm. 1.3] refined this
theorem by proving that any homogeneous polynomial relation between the values
$f_{1}(\alpha),\ldots,f_{n}(\alpha)$ with coefficients in
${\overline{\mathbb{Q}}}$ is a specialization of a homogeneous polynomial
relation between the functions $f_{1}(z),\ldots,f_{n}(z)$ with coefficients in
${\overline{\mathbb{Q}}}(z)$, again when $\alpha$ is not a singularity of $A$.
A less precise version of this theorem (but for $E$-functions in Siegel’s more
general sense; see Section 3.5 for details) had been obtained in 1996 by
Nesterenko and Shidlovskii [35], where $\alpha$ is simply assumed not to lie
in a certain finite set $S$, depending on the $f_{j}$’s but not specified in
their article. A fundamental consequence of their result is that a
transcendental $E$-function $f$ takes only finitely many algebraic values when
evaluated at algebraic points. To see this, one considers a non-trivial
minimal inhomogeneous differential equation with polynomial coefficients
$p+\sum_{j=0}^{\mu}p_{j}f^{(j)}=0$ satisfied by $f$ over
${\overline{\mathbb{Q}}}(z)$ and applies the Nesterenko-Shidlovskii theorem to
the functions $f_{1}:=1$, $f_{2}:=f,\ldots,f_{\mu}:=f^{(\mu-1)}$. They are
linearly independent over ${\overline{\mathbb{Q}}}(z)$, hence the numbers
$1,f(\alpha),\ldots,f^{(\mu-1)}(\alpha)$ are
${\overline{\mathbb{Q}}}$-linearly independent over ${\overline{\mathbb{Q}}}$
for all $\alpha\in{\overline{\mathbb{Q}}}\setminus S$; in particular
$f(\alpha)\notin{\overline{\mathbb{Q}}}$ for all such $\alpha$’s. For
$E$-functions in the strict sense, we now know thanks to Beukers [14] that if
$\alpha\in{\overline{\mathbb{Q}}}^{*}$ is not a root of the leading
coefficient $p_{\mu}$ above, then $f(\alpha)\notin{\overline{\mathbb{Q}}}$.

Algorithm 3 Algebraic values of $E$-functions over $\mathbb{Q}$

$\mathcal{L}=a_{r}(z)\partial_{z}^{r}+\dots+a_{0}(z)$;

Inputini: $f_{0}$ a truncated power series at precision $p_{0}\geq r$

specifying a unique solution $f\in\mathbb{Q}[[z]]$ of $\mathcal{L}(f)=0$.

Input _It is assumed that $f$ is an $E$-function._

Either “$f$ is a polynomial”,

or the finite set of all identities $f(\alpha)=\beta$ with algebraic $\alpha$
and $\beta$.

$\mathcal{L}_{\min}:=$MinimalRightFactor($\mathcal{L}$,ini)//  _Algorithm 1_

$\mathcal{L}_{\text{inhom}},g:=$MinimalInhomogeneousRightFactor($\mathcal{L}_{\min}$)//
_Algorithm 2_

if ${\operatorname{ord}}\mathcal{L}_{\text{inhom}}=0$ then return $f$ is a
polynomial

Define the polynomials $v_{0},\dots,v_{s+1}$ by
$\mathcal{L}_{\text{inhom}}(f)-g=v_{0}f^{(s)}-v_{1}f^{(s-1)}-\dots-
v_{s}f-v_{s+1}$

Form the companion matrix $M$ s.t.
$(0,f^{\prime},f^{\prime\prime},\dots,f^{(s)})^{\mathsf{T}}=M\cdot(1,f,f^{\prime},\dots,f^{(s-1)})^{\mathsf{T}}$

$\mathcal{R}:=\\{f(0)=f_{0}(0)\\}$

for $\mu\in{\mathbb{Q}}[z]$ irreducible factor of $v_{0}$ with $\mu(0)\neq 0$
do

Write $\alpha$ for a root of $\mu$

$\mathcal{M}:=$BeukersAlgo($M,\alpha$)

$(b_{1},\dots,b_{m}):=$basis of the left kernel of $\mathcal{M}(\alpha)$//
_Basis of algebraic relations at $\alpha$_

if there exists $(\beta,-1,0,\dots,0)$ in
$\operatorname{Span}(b_{1},\dots,b_{m})$ then

$\mathcal{R}:=\mathcal{R}\cup\\{f(\alpha)=\beta\\}$ return $\mathcal{R}$

Thus, in order to completely determine when an $E$-function takes a
transcendental value at a given non-zero algebraic point, one issue needs to
be dealt with: what happens for the (finite number of) algebraic numbers that
are roots of $p_{\mu}$ (in the same setting as above). This was done by
Adamczewski and Rivoal [5] by pushing further Beukers’ ideas from [14]. The
end result is an algorithm that takes as input an $E$-function $f$ and either
detects that $f$ is algebraic (in which case it is even a polynomial), or
computes the (finite) list of identities $f(\alpha)=\beta$ for algebraic
values $\alpha$ and $\beta$.

Algorithm 3 gives a version of the algorithm by Adamczewski and Rivoal that
benefits from the minimization algorithms of Section 2. It is stated here for
the $E$-functions over $\mathbb{Q}$ considered in this section (see the
comments in Sections 3.4 and 3.5 below for its extension to more general
settings). The algorithm relies on two results due to Beukers [14]:

  1. (1)

If $F=(f_{1},\dots,f_{n})^{\mathrm{T}}$ with $E$-functions $f_{i}$ is a
solution of $Y^{\prime}=AY$, the entries of $A$ belonging to
${\overline{\mathbb{Q}}}(z)$ and if $f_{1}(z),\dots,f_{n}(z)$ are linearly
independent over ${\overline{\mathbb{Q}}}(z)$, then for any non-zero $\alpha$
that is not a pole of $A$, the numbers $f_{1}(\alpha),\dots,f_{n}(\alpha)$ are
linearly independent over ${\overline{\mathbb{Q}}}$ [14, Corollary 1.4];

  2. (2)

Under the same assumptions, there exists a matrix $M$ with entries in
${\overline{\mathbb{Q}}}[z]$ such that $F=ME$, and $E$ is a vector of
$E$-functions solution to a system $Y^{\prime}=BY$ where $B$ does not have any
non-zero pole [14, Theorem 1.5].

Starting from a linear differential operator $\mathcal{L}$ and initial
conditions specifying an $E$-function $f$, the algorithm first computes a
minimal inhomogeneous equation of order $s$ for $f$. (This step also allows to
detect and discard a polynomial $f$.) By minimality of this equation,
$F=(1,f,f^{\prime},\dots,f^{(s-1)})$ is a vector of
${\overline{\mathbb{Q}}}(z)$-linearly independent $E$-functions solution to a
matrix deduced from the equation. Given a matrix $M$ as in the result (2)
above, it follows that the points $\alpha$ where
$(1,f(\alpha),\dots,f^{(s-1)}(\alpha))$ are linearly dependent over
${\overline{\mathbb{Q}}}$ are non-zero poles of $A$ where the left-kernel of
$M$ is not reduced to 0. The specific case of $f(\alpha)$ being algebraic
corresponds to the existence of a non-zero vector in that kernel whose first
two coordinates only are not zero. The remaining question is the computation
of these matrices $M$, which is described in Section 3.2.

Note that the first two steps of Algorithm 3, i.e. the calls to
MinimalRightFactor($\mathcal{L}$,ini) and
MinimalInhomogeneousRightFactor($\mathcal{L}_{\min}$), are not specific to
$E$-functions, and both return an output even when $f$ is not an $E$-function.
In that case, BeukersAlgo($M,\alpha$) terminates (by design) and it may even
output $\alpha$’s such that $f(\alpha)\in{\overline{\mathbb{Q}}}$; but it can
no longer be claimed that all such $\alpha$’s have been found.

###  3.2. Beukers’ algorithm and desingularization

Algorithm 3 concludes with a call to Algorithm BeukersAlgo($M,\alpha$)
described below. It is a clever desingularization process, which is different
from the one developed by Barkatou and Maddah [10], in that it does not rely
on Moser’s reduction [34, 9]. The end result is the following (Theorem 1.5 in
[14]).

######  Theorem 3.1.

Let $Y=(f_{1},\ldots,f_{n})^{\mathrm{T}}$ be a vector of
$\overline{{\mathbb{Q}}}(z)$-linearly independent $E$-functions satisfying
$Y^{\prime}=AY$, where $A$ is an $n\times n$ matrix with entries in
$\overline{{\mathbb{Q}}}(z)$. Then, there exists a vector of $E$-functions
$Z=(e_{1},\ldots,e_{n})^{\mathrm{T}}$ solution of $Z^{\prime}=BZ$ with $B$
having entries in $\overline{{\mathbb{Q}}}[z,1/z]$, and there exists a
polynomial matrix $M$ with entries in $\overline{{\mathbb{Q}}}[z]$ and
$\det(M)\neq 0$, such that
$(f_{1},\ldots,f_{n})^{\mathrm{T}}=M\cdot(e_{1},\ldots,e_{n})^{\mathrm{T}}$.

The key properties used in the proof of Theorem 3.1 are the statements (P1),
(P1’) and (P2) listed below. Note that for an $E$-function
$f\in{\mathbb{Q}}[[z]]$, we write $L^{\textrm{min}}_{f}$ for the monic linear
differential operator in ${\mathbb{Q}}(z)\langle\partial_{z}\rangle$ of
minimal order that cancels $f$.

(P1) For any $E$-function $f\in{\mathbb{Q}}[[z]]$, the finite non-zero
singularities of $L^{\textrm{min}}_{f}$ are apparent;

(P1’) If $A$ is an $n\times n$ matrix with entries in
$\overline{{\mathbb{Q}}}(z)$, and if $F$ is a vector of
$\overline{{\mathbb{Q}}}(z)$-linearly independent $E$-functions satisfying
$F^{\prime}=AF$, then the finite non-zero singularities of the system
$Y^{\prime}=AY$ are apparent;

(P2) If an $E$-function $f$ and $\alpha\in{\overline{\mathbb{Q}}}$ are such
that $f(\alpha)\in{\overline{\mathbb{Q}}}$, then
$({f(z)-f(\alpha)})/({z-\alpha})$ is an $E$-function.

Property (P1) is _André’s theorem_ [6, Cor. 4.4] and (P2) is an important
property of $E$-functions proved by Beukers [14, Prop. 4.1]. Property (P1’) is
a system version of André’s theorem, which is not, to our knowledge,
explicitly stated in the literature, although it is implicitly contained in
Beukers’ proof of his Theorem 1.5 in [14, p. 378]. For completeness, we detail
the proof of (P1’), which goes along the following lines.

_Proof of(P1’)_. Let $\mathcal{G}$ be the differential Galois group of (the
Picard-Vessiot field of) $Y^{\prime}=AY$. Let $V$ be the
${\overline{\mathbb{Q}}}$-vector space generated by the orbit
$\\{\sigma(F)\mid\sigma\in\mathcal{G}\\}$, where
$F=(f_{1},\ldots,f_{n})^{\mathrm{T}}$ is a vector of $E$-functions satisfying
$F^{\prime}=AF$, linearly independent over ${\overline{\mathbb{Q}}}(z)$. The
conclusion of (P1’) clearly follows by combining the following two steps.

Step 1. The dimension of $V$ over ${\overline{\mathbb{Q}}}$ is equal to $n$,
hence one can extract from $V$ a fundamental matrix of solutions
$\mathcal{F}$, whose first column is $F$.

Step 2. All the entries of $\mathcal{F}$ are holomorphic at all non-zero
points $\alpha\in\mathbb{C}\setminus\\{0\\}$.

_Proof of Step 2_: Let $L_{i}:=L^{\textrm{min}}_{f_{i}}$. Since elements of
$\mathcal{G}$ commute with differentiation, all $\sigma(f_{i})$ are solutions
of $L_{i}$ for all $\sigma\in\mathcal{G}$. By André’s theorem (P1),
$\sigma(f_{i})$ has no true singularity in $\mathbb{C}\setminus\\{0\\}$. Hence
$\mathcal{F}$ is holomorphic at any $\alpha\in\mathbb{C}\setminus\\{0\\}$.

_Proof of Step 1_: If $A$ is a companion matrix, then the shape of the system
$Y^{\prime}=AY$ implies that $F$ is of the form
$F=(f,f^{\prime},\ldots,f^{(n-1)})^{\mathrm{T}}$, where $f$ is the
$E$-function $f=f_{1}$. The linear independence assumption implies that
$L^{\textrm{min}}_{f}$ has order $n$. By [49, Corollary 1.38] (see also [13,
p. 190, Proposition 3]), the dimension of the ${\overline{\mathbb{Q}}}$-vector
space $\tilde{V}$ generated by the orbit
$\\{\sigma(f)\mid\sigma\in\mathcal{G}\\}$ is equal to $n$. On the one hand,
this dimension is upper bounded by the dimension of $V$, since any linear
relation among the entries of $V$ yields a linear relation among the elements
of the set $\\{\sigma(f)\mid\sigma\in\mathcal{G}\\}$. On the other hand, $V$
is included in the solution space of $Y^{\prime}=AY$, hence it has dimension
at most $n$. Therefore, $\dim_{{\overline{\mathbb{Q}}}}(V)=n$, and the
assertion is proved in the companion case.

Now, if $A$ a general matrix, by the cyclic vector lemma (see e.g. [18, Thm
3.11], or [49, Proposition 2.9]) the system $Y^{\prime}=AY$ is “gauge
equivalent” to $Z^{\prime}=CZ$, where $C$ is a companion matrix with entries
in ${\overline{\mathbb{Q}}}(z)$. This means that there exists an invertible
matrix $P$ with entries in ${\overline{\mathbb{Q}}}(z)$ such that $Z:=P\cdot
Y$ satisfies $Z^{\prime}=P[A]\cdot Z$, where $P[A]:=(PA+P^{\prime})P^{-1}$ is
equal to a companion matrix $C$. Moreover, by [19, §6] (see also [49, Lemma
2.10]), the entries of the matrix $P$ can be chosen to be polynomials in
${\overline{\mathbb{Q}}}[z]$, of degree at most $n-1$. By construction,
$G:=P\cdot F$ satisfies $G^{\prime}=C\cdot G$. Hence, the vector $G$ is
necessarily of the form $G:=(g,g^{\prime},\ldots,g^{(n-1)})^{\mathrm{T}}$,
where $g$ is a ${\overline{\mathbb{Q}}}[z]$-linear combination of the
$E$-functions $f_{i}$. In particular, $g$ is itself an $E$-function. Moreover,
$L^{\textrm{min}}_{g}$ has order $n$: indeed, any
${\overline{\mathbb{Q}}}[z]$-linear combination $0=\mathbf{v}\cdot G$ between
$g,g^{\prime},\ldots,g^{(n-1)}$ yields a ${\overline{\mathbb{Q}}}[z]$-linear
combination $0=(\mathbf{v}P)\cdot F$ between the entries of $F$; since these
are assumed linearly independent over ${\overline{\mathbb{Q}}}(z)$, and since
$P$ is invertible, $\mathbf{v}$ is necessarily zero. We are now in position to
apply the companion case. Since gauge equivalent systems have the same
differential Galois group [45, p. 13], the new companion system
$Z^{\prime}=CZ$ has differential Galois group $\mathcal{G}$. By applying the
companion case, we deduce that $\dim_{{\overline{\mathbb{Q}}}}(V_{C})=n$,
where $V_{C}$ is the ${\overline{\mathbb{Q}}}$-vector space generated by the
orbit $\\{\sigma(G)\,|\,\sigma\in\mathcal{G}\\}$. It remains to show that
$\dim_{{\overline{\mathbb{Q}}}}(V)=\dim_{{\overline{\mathbb{Q}}}}(V_{C})$.
Choose $\sigma_{1},\ldots,\sigma_{n}$ in $\mathcal{G}$ such that
$\sigma_{1}(G),\ldots,\sigma_{n}(G)$ are linearly independent over
${\overline{\mathbb{Q}}}$. Then, $\sigma_{1}(F),\ldots,\sigma_{n}(F)$ are also
linearly independent over ${\overline{\mathbb{Q}}}$, because of the relation
$G=P\cdot F$ and the fact that all elements in $\mathcal{G}$ leave $P$
invariant. It follows that $V$ has dimension at least $n$; since $V$ is
included in the solution space of $Y^{\prime}=AY$, it also has dimension at
most $n$, therefore $\dim_{{\overline{\mathbb{Q}}}}(V)=n$, which concludes the
proof. $\square$

We now prove Theorem 3.1. Our proof is inspired by that of Beukers in [14].
The main difference is that our proof does not depend on a specific
desingularization procedure for linear differential systems.

###### Proof of Theorem 3.1.

We make use of a _desingularization lemma_ [10, Theorem 2]: there exists a
polynomial matrix $M$ with entries in $\overline{{\mathbb{Q}}}[z]$ and with
$\det(M)\neq 0$ such that the finite poles of $B=M[A]:=M^{-1}(AM-M^{\prime})$
are exactly the true (i.e, non-apparent) singularities of $Y^{\prime}=AY$ and
such that $\det(M)$ is a non-zero polynomial whose roots are among the
apparent singularities of $Y^{\prime}=AY$. (See also Proposition 3.2 below.)

In our case, by Property (P1’) above, the entries of $B$ are in
$\overline{{\mathbb{Q}}}[z,1/z]$.

Define
$Z=(e_{1},\ldots,e_{n})^{\mathrm{T}}:=M^{-1}\cdot(f_{1},\ldots,f_{n})^{\mathrm{T}}$,
so that
$(f_{1},\ldots,f_{n})^{\mathrm{T}}=M\cdot(e_{1},\ldots,e_{n})^{\mathrm{T}}$. A
simple computation shows that $Z^{\prime}=BZ$. It remains to prove that all
the $e_{i}$’s are $E$-functions. The proof relies on Property (P2) above.

By definition, each $e_{i}$ is equal to
$\frac{1}{\det(M)}\cdot\sum_{j=1}^{n}p_{i,j}f_{j}$ for polynomials $p_{i,j}$
in $\overline{{\mathbb{Q}}}[z]$. Since $B$ has no non-zero pole, each $e_{i}$
is holomorphic at every apparent singularity $\rho\neq 0$ of $Y^{\prime}=AY$.
Therefore, $\sum_{j=1}^{n}p_{i,j}f_{j}$ is an $E$-function which vanishes at
any root of $\det(M)$ at an order at least equal to the multiplicity of that
root in $\det(M)$. By repeated application of Property (P2), it follows that
$e_{i}$ is an $E$-function. ∎

Beukers’ proof [14, p. 378] of Theorem 3.1 actually contains a general
effective desingularization process, which deserves to be stated independently
of the context of $E$-functions. It is given in Algorithm 4, whose properties
are summarized in the following.

######  Proposition 3.2.

Let $A$ be an $n\times n$ matrix with entries in ${\mathbb{Q}}(z)$ and let
$\alpha\in\overline{\mathbb{Q}}$ be such that a fundamental solution
$\mathcal{Y}$ of $Y^{\prime}=AY$ is holomorphic at $\alpha$. Then Algorithm 4
computes a matrix of polynomials $B\in(\mathbb{Q}(\alpha)[z])^{n\times n}$
such that $\mathcal{Y}=B\mathcal{Z}$ with $\mathcal{Z}$ a fundamental solution
of $Z^{\prime}=CZ$ also holomorphic at $\alpha$ and
$C\in(\mathbb{Q}(\alpha)(z))^{n\times n}$ only has poles where $A$ does,
except at $\alpha$, where it is holomorphic.

###### Proof.

We reproduce Beukers’ proof, with more details.

By hypothesis, the determinant $W=\det\mathcal{Y}$ is holomorphic in a
neighborhood of $z=\alpha$. If $W(\alpha)\neq 0$, then $\mathcal{Y}^{-1}$ is
holomorphic in the neighborhood of $z=\alpha$ and therefore so is
$A=\mathcal{Y}^{\prime}\mathcal{Y}^{-1}$. In that case, $C=A$ and
$B={\operatorname{Id}}_{n}$ gives the result.

Otherwise, as $W\neq 0$, there exists $r\in\mathbb{N}_{>0}$ such that
$W(z)\sim c(z-\alpha)^{r}$ for $z\rightarrow\alpha$ with $c\neq 0$. Since $W$
satisfies $W^{\prime}=\operatorname{Trace}(A)W$, it follows that $r$ is the
residue of $\operatorname{Trace}(A)$ at $z=\alpha$. Starting with
$B={\operatorname{Id}}_{n}$, $C=A$ and $\mathcal{Z}=\mathcal{Y}$, the
algorithm repeats at most $r$ times an operation that updates $C$ and $B$ so
that $\mathcal{Y}=B\mathcal{Z}$ and

  * —

$B$ is a matrix of polynomials;

  * —

$\mathcal{Z}$ is holomorphic at $\alpha$;

  * —

$C:=B^{-1}(AB-B^{\prime})$ has no pole outside those of $A$;

  * —

$\mathcal{Z}$ is a fundamental solution of $Z^{\prime}=CZ$;

  * —

$\operatorname{val}_{z=\alpha}\det\mathcal{Z}=\operatorname{val}_{z=\alpha}\det\mathcal{Y}-1$.

By composing these steps, it is sufficient to prove that one iteration of the
loop has these properties.

Each step is centered around the definition of a matrix $M$ as follows. Let
$k>0$ be the order of the pole of the matrix $C$ at $\alpha$, let $i$ be the
index of the first row of $C$ with a pole of order $k$, let $v$ be the
constant vector of coefficients of $(z-\alpha)^{-k}$ in that row, $D$ be the
diagonal matrix $\operatorname{diag}(1,\dots,1,z-\alpha,1,\dots,1)$ with
$z-\alpha$ in the $i$th position and $M$ be an invertible constant matrix with
$v$ in its $i$th row. Then $\tilde{B}:=BM^{-1}D$ possesses the desired
properties:

  * —

$\tilde{B}$ is the product of matrices of polynomials;

  * —

$\tilde{\mathcal{Z}}:=D^{-1}M\mathcal{Z}$ is holomorphic at $\alpha$: since
both $\mathcal{Z}$ and $\mathcal{Z}^{\prime}$ are holomorphic at $\alpha$, the
product $v\mathcal{Z}$ is 0 at $\alpha$, making the $i$th row of
$M\mathcal{Z}$ a multiple of $(z-\alpha)$;

  * —

$\tilde{C}:=D^{-1}M(CM^{-1}D-M^{-1}D^{\prime})$ has no pole outside those of
$C$;

  * —

$\tilde{\mathcal{Z}}$ is a fundamental solution of $Z^{\prime}=\tilde{C}Z$;

  * —

$\det\tilde{\mathcal{Z}}=\det M\det\mathcal{Z}/(z-\alpha).$∎

Algorithm 4 Removal of Singularities (BeukersAlgo($M,\alpha$))

$A$: matrix in ${\mathbb{Q}}(z)^{n\times n}$;

Input$\alpha$: root of the denominator of an entry in $A$

$B$: matrix in $({\mathbb{Q}}(\alpha)[z])^{n\times n}$ satisfying Proposition
3.2.

$r:=\operatorname{Res}_{z=\alpha}\operatorname{Trace}(A)$//  _Residue of the
trace_

if $r\not\in\mathbb{N}_{\geq 0}$ then error singularity cannot be removed

$C:=A$; $B:={\operatorname{Id}}_{n}$

for $m=1,\dots,r$ do

$k:=$order of the pole of $C$ at $\alpha$

if $k=0$ then break

$i:=$index of the first row of $C$ with a pole of order $k$ at $\alpha$

$v:=$vector of coefficients of $(z-\alpha)^{-k}$ in row $i$ of $C$

$M:=$an invertible constant matrix with $v$ its $i$th row; //  _Complete $v$
into a basis_

$D:=\operatorname{diag}(1,\dots,1,z-\alpha,1,\dots,1)$ with $z-\alpha$ in the
$i$th position

$B:=BM^{-1}D$

$C:=D^{-1}MCM^{-1}D-D^{-1}D^{\prime}$

return $B$

###  3.3. Effective decomposition of $E$-functions

For an $E$-function $f\in{\mathbb{Q}}[[z]]$ (or more generally in
${\overline{\mathbb{Q}}}[[z]]$), we call _exceptional values_ those (finitely
many) non-zero algebraic numbers $\alpha$ such that
$f(\alpha)\in{\overline{\mathbb{Q}}}$. The set of exceptional values of $f$ is
denoted by $\textrm{Exc}(f)$. We call the $E$-function $f$ _purely
transcendental_ if it has no exceptional values, i.e. if
$\textrm{Exc}(f)=\emptyset$.

This subsection deals with the fact that every $E$-function is equal to the
sum of a polynomial and of a polynomial multiple of a purely transcendental
$E$-function. The existence of such a decomposition was proved in [41] for
general $E$-functions, not necessarily with coefficients in ${\mathbb{Q}}$
(see Section 3.4). Such a decomposition is not unique: indeed, if we have
$f=p+qg$ where $p,q$ are in ${\overline{\mathbb{Q}}}[z]$ and $g$ is a purely
transcendental $E$-function, then for any $u\in{\overline{\mathbb{Q}}}[z]$,
the identity $f=p-qu+q(g+u)$ is another admissible decomposition because $g+u$
is still a purely transcendental $E$-function. In particular, the Taylor
coefficients of $f$ and $g$ may lie in two different number fields. However,
we have the following:

######  Proposition 3.3.

Every transcendental $E$-function (with coefficients in
${\overline{\mathbb{Q}}}$) can be written in a unique way as $f=p+qg$ with
$p,q\in{\overline{\mathbb{Q}}}[z]$, $q$ monic and $q(0)\neq 0$,
$\deg(p)<\deg(q)$ and $g$ a purely transcendental $E$-function.

We shall call this decomposition the canonical decomposition of $f$.

###### Proof.

Consider two decompositions $p+qg=\widetilde{p}+\widetilde{q}\widetilde{g}$ of
an $E$-function $f$. Since $f$ is transcendental, $q,g,\widetilde{q}$ and
$\widetilde{g}$ are not identically zero. We can and shall assume that
$q(0)\widetilde{q}(0)\neq 0$, because we can always multiply $g$ and
$\widetilde{g}$ by a suitable power of $z$ and keep the “purely
transcendental” property. This being done, obviously $q$ and $\widetilde{q}$
share the same set of roots, which is $\textrm{Exc}(f)$. Moreover $g$ and
$\widetilde{g}$ being purely transcendental, we claim that any root of $q$ and
$\widetilde{q}$ has the same multiplicity in $q$ and in $\widetilde{q}$, so
that $q$ and $\widetilde{q}$ are equal up to a non-zero constant factor. To
prove the claim, let $\rho$ be a root of $q$ of multiplicity $m$ and of
multiplicity $\widetilde{m}$ for $\widetilde{q}$: if $m<\widetilde{m}$ then
differentiating $m$ times both sides of
$p+qg=\widetilde{p}+\widetilde{q}\widetilde{g}$ and evaluating at $z=\rho$, we
obtain that $g(\rho)\in{\overline{\mathbb{Q}}}$ which is not possible because
$\rho\neq 0$, hence by symmetry of the situation we have $m=\widetilde{m}$.
Consequently choosing $q$ monic makes it unique.

Given two decompositions $p+qg=\widetilde{p}+q\widetilde{g}$ of $f$, we have
$g-\widetilde{g}\in{\overline{\mathbb{Q}}}[z]$. Indeed, the $E$-function
$g-\widetilde{g}=(\widetilde{p}-p)/q$ is a rational function hence a
polynomial. Therefore, by Euclidean division of $p$ by $q$, this property
implies that it is always possible to impose $\deg(p)<\deg(q)$.

Now, if in the decompositions $p+qg=\widetilde{p}+q\widetilde{g}$, we have
$\deg(p)<\deg(q)$ and $\deg(\widetilde{p})<\deg(q)$, then necessarily
$p=\widetilde{p}$ and $g=\widetilde{g}$: indeed, we have
$g-\widetilde{g}=u\in{\overline{\mathbb{Q}}}[z]$, so that
$p-\widetilde{p}=-uq$ which by consideration of the degree on both sides
forces $u=0$. ∎

Our main contribution is to prove that when $f$ has coefficients in
$\mathbb{Q}$, then we can find polynomials and a purely transcendental
$E$-function involved in the canonical decomposition that also have
coefficients in $\mathbb{Q}$, and moreover that one can compute this
decomposition algorithmically. More precisely, we prove:

######  Theorem 3.4.

Any transcendental $E$-function $f\in{\mathbb{Q}}[[z]]$ admits a canonical
decomposition $f=p+qg$, where $p$ and $q$ are polynomials in ${\mathbb{Q}}[z]$
and $g\in{\mathbb{Q}}[[z]]$ is a purely transcendental $E$-function. Moreover,
if $f$ is given by a linear differential equation together with sufficiently
many initial terms, then one can effectively determine $p$ and $q$.

Of course, once $p$ and $q$ are determined, $g$ is determined by a linear
differential equation together with sufficiently many initial terms, simply
because $g=(f-p)/q$. Before proceeding to the proof of Theorem 3.4, we state a
very useful fact.

######  Proposition 3.5.

Let $f\in{\mathbb{Q}}[[z]]$ be an $E$-function and let
$\alpha\in\mathrm{Exc}(f)$. Then,

  * $(i)$

$f(\alpha)\in\mathbb{Q}(\alpha)$;

  * $(ii)$

all Galoisian conjugates of $\alpha$ belong to $\mathrm{Exc}(f)$;

  * $(iii)$

for any Galoisian conjugate $\alpha^{\prime}$ of $\alpha$, the value
$f(\alpha^{\prime})$ is a Galoisian conjugate of $f(\alpha)$.

###### Proof.

The three statements are consequences of Algorithm 3. For any given root
$\alpha$ of any irreducible factor $\mu\in\mathbb{Q}[z]$ of $v_{0}$, the
algorithm determines if there exists a vector in
${\overline{\mathbb{Q}}}^{s+1}$ of the form $(\beta,-1,0,\ldots,0)$ which is
in the left kernel of the matrix $\mathcal{M}(\alpha)$, whose entries are in
$\mathbb{Q}(\alpha)$. The existence of this vector is equivalent to
$f(\alpha)=\beta\in{\overline{\mathbb{Q}}}$. When it exists, this proves that
$f(\alpha)=\beta\in\mathbb{Q}(\alpha)$, a fact proved in [22, 23] in the
general case (with $\mathbb{Q}(\alpha)$ replaced by $\mathbb{K}(\alpha)$ when
$f$ has coefficients in a number field $\mathbb{K}$). Now, such a vector
exists if and only there exists a vector of the same form in the left kernel
of the matrix $\mathcal{M}(\alpha^{\prime})$, where $\alpha^{\prime}$ is any
Galoisian conjugate of $\alpha$ (i.e. any other root of $\mu$ in this case).
It follows that for any conjugate $\alpha^{\prime}$ of $\alpha$,
$f(\alpha^{\prime})$ is a conjugate of $f(\alpha)$. ∎

###### Proof of Theorem 3.4.

The set $\textrm{Exc}(f)$ can be computed using Algorithm 3. If
$\textrm{Exc}(f)=\emptyset$, then the canonical decomposition of $f$ is
$f=p+qg$ with $p=0$, $q=1$ and $g=f$. From now on, we assume that
$\textrm{Exc}(f)=\\{\alpha_{1},\ldots,\alpha_{k}\\}\neq\emptyset$; by
Proposition 3.5, $\textrm{Exc}(f)$ can be partitioned into blocks of Galois
conjugated values.

Let us first assume that there is only one such block, i.e. that
$\\{\alpha_{1},\ldots,\alpha_{k}\\}$ is the set of roots of a monic
irreducible polynomial $E\in{\mathbb{Q}}[z]$. For any $m\geq 0$, we consider
the $E$-adic expansion of $f$ to order $m$:

$f=p_{0}+p_{1}E+\cdots+p_{m}E^{m}+E^{m+1}g_{m},$ (13)

with $p_{0}(z),\ldots,p_{m}(z)\in\mathbb{C}[z]$ each of degree less than
$k=\deg(E)$, and $g_{m}\in\mathbb{C}[[z]]$.

We will prove the following claims:

Claim 1. $p_{0},\ldots,p_{m}\in{\mathbb{Q}}[z]$ and
$g_{m}\in{\mathbb{Q}}[[z]]$.

Claim 2. $g_{m}$ is an $E$-function.

Claim 3. There exists an $m\geq 0$ such that $g_{m}$ is purely transcendental.

From these claims, the proof of the first part of the theorem follows by
taking $p:=p_{0}+p_{1}E+\cdots+p_{m}E^{m}$, $q:=E^{m+1}$ and $g:=g_{m}$.

Proof of Claim 1. It is enough to prove it for $m=0$, and then iterate. We
have $f=p_{0}+Eg_{0}$ with $p_{0}\in\mathbb{C}[z]$ of degree less than $k$ and
$g_{0}\in\mathbb{C}[[z]]$ and we need to prove that the coefficients of
$p_{0}$ and $g_{0}$ are actually in ${\mathbb{Q}}$. First, we observe that
$p_{0}(z)=c_{0}+c_{1}z+\cdots+c_{k-1}z^{k-1}$ is the unique polynomial in
$\mathbb{C}[z]$ such that $p_{0}(\alpha_{i})=f(\alpha_{i})$ for $1\leq i\leq
k$. In matrix terms, this rewrites as

$\begin{pmatrix}1&\alpha_{1}&\ldots&\alpha_{1}^{k-1}\\\ \vdots&&&\vdots\\\
1&\alpha_{k}&\ldots&\alpha_{k}^{k-1}\\\
\end{pmatrix}\cdot\begin{pmatrix}c_{0}\\\ \vdots\\\ c_{k-1}\\\
\end{pmatrix}=\begin{pmatrix}f(\alpha_{1})\\\ \vdots\\\ f(\alpha_{k-1})\\\
\end{pmatrix}$

and by multiplying this equality on the left by the transpose of the
Vandermonde matrix, we get the equivalent identity

$\begin{pmatrix}k&\sum_{i}\alpha_{i}&\ldots&\sum_{i}\alpha_{i}^{k-1}\\\
\sum_{i}\alpha_{i}&\sum_{i}\alpha_{i}^{2}&\ldots&\sum_{i}\alpha_{i}^{k}\\\
\vdots&&&\vdots\\\
\sum_{i}\alpha_{i}^{k-1}&\sum_{i}\alpha_{i}^{k}&\ldots&\sum_{i}\alpha_{i}^{2k-2}\\\
\end{pmatrix}\cdot\begin{pmatrix}c_{0}\\\ c_{1}\\\ \vdots\\\ c_{k-1}\\\
\end{pmatrix}=\begin{pmatrix}\sum_{i}f(\alpha_{i})\\\
\sum_{i}\alpha_{i}f(\alpha_{i})\\\ \vdots\\\
\sum_{i}\alpha_{i}^{k-1}f(\alpha_{i})\\\ \end{pmatrix}.$ (14)

Now, the matrix on left-hand side of Eq. 14 is invertible and with
coefficients in ${\mathbb{Q}}$, since it contains the power sums of the roots
of the polynomial $E\in{\mathbb{Q}}[z]$. On the other hand, for any
$E$-function $g\in{\mathbb{Q}}[[z]]$, we have that
$\sum_{i}g(\alpha_{i})\in{\mathbb{Q}}$, by Proposition 3.5. Applying this to
the $E$-functions $f(z)$, $zf(z),\ldots,z^{k-1}f(z)$, we deduce that the
right-hand side of Eq. 14 is a vector in ${\mathbb{Q}}^{k}$. This implies that
the $c_{i}$’s are all in ${\mathbb{Q}}$, hence $p_{0}\in{\mathbb{Q}}[z]$. From
there it directly follows that $g_{0}\in{\mathbb{Q}}[[z]]$.

Proof of Claim 2. It is again enough to prove the claim for $m=0$, and then
iterate. Indeed, from Eq. 13 it follows that the $E$-adic expansion of
$g_{m-1}$ to order 1 is $g_{m-1}=p_{m}+Eg_{m}$, and since $f=p_{0}+Eg_{0}$ is
an $E$-function one deduces iteratively that $g_{0},g_{1},\ldots,g_{m}$ are
$E$-functions.

It remains to prove that if we have $f=p+Eg$ with
$E\in{\overline{\mathbb{Q}}}[z]$ and $p\in{\overline{\mathbb{Q}}}[z]$ of
degree less than $k=\deg(E)$ and $g\in{\overline{\mathbb{Q}}}[[z]]$, then $g$
is an $E$-function. This is done by induction on $k\geq 1$. For $k=1$, this is
precisely Property (P2). Assume the property is proved for any $E$ of degree
$k-1\geq 1$ and any $p$ of degree less than $k-1$. Assume we have $f=p+Eg$
with $E\in{\overline{\mathbb{Q}}}[z]$ of degree $k$,
$p\in{\overline{\mathbb{Q}}}[z]$ of degree less than $k$ and
$g\in{\overline{\mathbb{Q}}}[[z]]$. Let $\beta$ be one of the roots of $E$ and
write $p(z)=\sum_{j=0}^{k-1}p_{j}(z-\beta)^{j}$. Then $f(\beta)=p_{0}$ and by
Property (P2),

$\frac{f(z)-f(\beta)}{z-\beta}=\sum_{j=0}^{k-2}p_{j+1}(z-\beta)^{j}+\frac{E(z)}{z-\beta}g(z)$

is an $E$-function. Since $E(z)/(z-\beta)$ is a polynomial of degree $k-1$ and
$\sum_{j=0}^{k-2}p_{j+1}(z-\beta)^{j}$ is of degree less that $k-1$, we deduce
that $g$ is an $E$-function by the induction hypothesis.

Proof of Claim 3. From Eq. 13 it follows that the only exceptional values of
the $g_{m}$’s are necessarily contained in the set
$\textrm{Exc}(f)=\\{\alpha_{1},\ldots,\alpha_{k}\\}$.

We will show that there exists an $m$ such that $g_{m}$ does not have any of
the $\alpha_{j}$’s as an exceptional value, and therefore $g_{m}$ is purely
transcendental. By Proposition 3.5, this is equivalent to proving that the
$g_{m}$’s cannot all share $\alpha:=\alpha_{1}$ as an exceptional value.

Setting $g_{-1}=f$, it follows from Eq. 13 (with $m$ replaced by $m-1$, and
then by differentiating $m$ times) that $1,f^{(m)}(\alpha)$ and
$g_{m-1}(\alpha)$ are linearly dependent over ${\overline{\mathbb{Q}}}$ for
all $m\geq 0$ by a relation of the form
$f^{(m)}(\alpha)=u_{m}+v_{m}g_{m-1}(\alpha)$ with
$u_{m},v_{m}\in{\overline{\mathbb{Q}}}$ and $v_{m}\neq 0$. Hence,

$\text{trdeg}_{{\overline{\mathbb{Q}}}}\big{(}f(\alpha),\ldots,f^{(m)}(\alpha)\big{)}=\text{trdeg}_{{\overline{\mathbb{Q}}}}\big{(}g_{-1}(\alpha),\ldots,g_{m-1}(\alpha)\big{)}\quad\text{for
all}\;m\geq 0.$

By contradiction, let us now assume that
$g_{m}(\alpha)\in{\overline{\mathbb{Q}}}$ for all $m\geq-1$. Then we have
$\text{trdeg}_{{\overline{\mathbb{Q}}}}(f(\alpha),\ldots,f^{(m)}(\alpha))=0$
for all $m\geq 0$. Now by Property (P1) it follows that $f$ satisfies an ODE,
of some order $\mu\geq 1$, having only $z=0$ as finite singularity. By
considering the corresponding companion system $Y^{\prime}=AY$ where $f$ is
the first element of the column vector $Y$, the matrix $A$ has Laurent
polynomial entries in $z$, hence the Siegel-Shidlovskii theorem ensures that

$0=\text{trdeg}_{{\overline{\mathbb{Q}}}}\big{(}f(\alpha),\ldots,f^{(\mu-1)}(\alpha)\big{)}=\text{trdeg}_{{\overline{\mathbb{Q}}}(z)}\big{(}f(z),\ldots,f^{(\mu-1)}(z)\big{)}\geq
1,$

a contradiction.

On the effective side, note that one can compute the $E$-adic expansion (13)
of $f$ to any order $m$, for instance using linear algebra. Then, to compute
the needed decomposition, one may, for increasing values $m=0,1,\ldots$,
compute a linear differential equation for $g_{m}$ as in (13) together with
sufficiently many initial terms, and test using Algorithm 3 whether
$\mathrm{Exc}(g_{m})$ is empty or not. This procedure will eventually
terminate.

We now treat the general case, where $\mathrm{Exc}(f)$ contains several blocks
$B_{1},\ldots,B_{p}$, each block containing conjugated exceptional values.
Denote by $E_{j}(z)$ the minimal polynomial $\prod_{\alpha\in
B_{j}}(z-\alpha)$ of the elements in $B_{j}$. By the reasoning used in the
case of a single block, one first finds a decomposition $f=p_{1}+q_{1}g_{1}$
with $p_{1},q_{1}$ in ${\mathbb{Q}}[z]$ and $g_{1}\in{\mathbb{Q}}[[z]]$ an
$E$-function such that $\mathrm{Exc}(g_{1})=\textrm{Exc}(f)\setminus B_{1}$.
Then, one applies the same to the $E$-function $g_{1}$, and writes it as
$g_{1}=p_{2}+q_{2}g_{2}$, and thus $f=(p_{1}+q_{1}p_{2})+(q_{1}q_{2})g_{2}$,
with $p_{2},q_{2}$ in ${\mathbb{Q}}[z]$ and $g_{2}\in{\mathbb{Q}}[[z]]$ an
$E$-function such that $\mathrm{Exc}(g_{2})=\textrm{Exc}(f)\setminus(B_{1}\cup
B_{2})$. Continuing the same way $p$ times, we end up with a decomposition
$f=p+qg$, with $p,q$ in ${\mathbb{Q}}[z]$ and $g\in{\mathbb{Q}}[[z]]$ an
$E$-function such that
$\mathrm{Exc}(g)=\textrm{Exc}(f)\setminus(B_{1}\cup\cdots\cup
B_{p})=\emptyset$. Moreover, by construction we have that $q$ monic, $q(0)\neq
0$ and $\deg(p)<\deg(q)$. This concludes the proof. ∎

###  3.4. $E$-functions with coefficients in a number field

In general, an $E$-function is a power series

$f(z):=\sum_{n=0}^{\infty}\frac{a_{n}}{n!}z^{n}\quad{\text{in}}\quad{\overline{\mathbb{Q}}}[[z]]$

such that

  1. $(i)$

$f(z)$ satisfies a homogeneous linear differential equation with coefficients
in ${\overline{\mathbb{Q}}}(z)$;

there exists $C>0$ such that

  1. $(ii)$

for any $\sigma\in\textup{Gal}({\overline{\mathbb{Q}}}/\mathbb{Q})$ and any
$n\geq 0$, $|\sigma(a_{n})|\leq C^{n+1}$;

  2. $(iii)$

for any $n\geq 0$, there exists $d_{n}\in\mathbb{N}\setminus\\{0\\}$ such that
$d_{n}\leq C^{n+1}$ and $d_{n}a_{m}\in\mathcal{O}_{{\overline{\mathbb{Q}}}}$
for all $0\leq m\leq n$.

In particular $(ii)$ with $\sigma=\textup{id}$ implies that $f(z)$ is an
entire function. Moreover, $(i)$ implies that the coefficients $a_{n}$ all
live in a certain number field, so that there are only finitely many Galoisian
conjugates to consider in $(ii)$; if $a_{n}\in\mathbb{Q}$, this definition
reduces to that of Section 3.1.

The Adamczewski-Rivoal algorithm applies to these more general situations. The
version stated in Section 3.1 also applies. Indeed, all the tools it uses work
more generally. This is obviously the case for Beukers’ desingularization, it
is also the case for the algorithms used by minimization: greatest common
right divisors, Hermite-Padé approximants, series solutions and the
computation of bounds on the degrees of factors (see [17]).

###  3.5. Siegel’s original definition

$E$-functions with algebraic coefficients have been first defined by Siegel
[44] in 1929 in a more general way: in $(ii)$ and $(iii)$ above, the upper
bounds $(\cdots)\leq C^{n+1}$ for all $n\geq 0$ are replaced by: for all
$\varepsilon>0$, there exists $N(\varepsilon)$ such that $(\cdots)\leq
n!^{\varepsilon}$ for all $n\geq N(\varepsilon)$. $E$-functions considered
above are sometimes denoted $E^{*}$-functions or called “$E$-functions in the
strict sense”: since André’s work [6, 7], it has become standard (though
improper) to call them simply “$E$-functions” as well.

The Siegel-Shidlovskii and Nesterenko-Shidlovskii theorems both hold in that
setting. The latter was refined by Beukers for $E$-functions in the strict
sense only. Then, André generalized Beukers’ lifting theorem to $E$-functions
in Siegel’s sense by a completely different method [8]; another proof was
later given by Lepetit [30] by a (non-trivial) adaptation of Beukers’ original
method.

We note here that Lepetit [30] also generalized the Adamczewski-Rivoal
algorithm to the case of $E$-functions in Siegel’s original sense: he showed
that all the steps in this algorithm work exactly the same mutatis mutandis,
so that in fact our more efficient algorithm described here applies as well if
the input is an $E$-function in Siegel’s sense with rational coefficients.
Moreover, the decomposition $f=p+qg$ studied in Section 3.3, holds in Siegel’s
sense, in particular Theorem 3.4. However, it is conjectured that the classes
of $E$-functions in Siegel’s sense and of $E$-functions in the strict sense
are the same (see [6, p. 715]). This implies that the distinction is in
practice illusory because all known examples of $E$-functions satisfy all the
conditions to be $E$-functions in the strict sense.

##  4\. Examples and implementation

###  4.1. Minimization is simpler than factorization

The following simple example illustrates the difference between minimization
and factorization. Take

$A=z^{2}\partial_{z}+3,\quad B=(z-3)\partial_{z}+4z^{5}$

and their product

$C=AB=z^{2}(z-3)\partial_{z}^{2}+(4z^{7}+z^{2}+3z-9)\partial_{z}+4z^{5}(5z+3).$

The computation of a bound on the degree of the coefficients of a factor of
order 1 of $C$ gives 970. This leads to a large computation when trying to
factor $C$ without further information. With the extra knowledge that we are
looking for a solution of $C$ with initial condition $y(0)=1$, we easily
compute the first 20 coefficients of the unique series solution $S$ of $C$
with $S(0)=1$ and then compute an approximant basis (a Hermite-Padé
approximant) of $(S,S^{\prime})$ at order 20. This recovers the operator $B$.
It is easily checked that $B$ is a right divisor of $C$ by Euclidean division.
It follows that all solutions of $B(y)=0$ such that $y(0)=1$ are solutions of
$C(y)=0$ and by uniqueness this proves that $B(y)=0$ is the minimal
homogeneous differential equation for $S$. This takes less than a second with
our implementation (111Minimization is available as part of the gfun package
available at https://perso.ens-lyon.fr/bruno.salvy/software/the-gfun-package/;
the code for exceptional algebraic values of $E$-functions can also be
downloaded from that page, together with a worksheet of examples.). On this
example, Maple’s factorization routine DEtools[DFactor] has to be killed after
running for 1 hour, a further indication that factorization is more complex
than minimization.

###  4.2. The Lorch-Muldoon example

In a special case of a result due to Lorch and Muldoon [32], the starting
point is the following equation satisfied by the fourth derivative of Bessel’s
$J_{0}$ function:

$(z^{5}-10z^{3}+45z)y(z)+(z^{4}-18z^{2}+45)y^{\prime}(z)+(z^{5}-6z^{3}+9z)y^{\prime\prime}(z)=0,$

with initial conditions $y(0)=3/8,y^{\prime}(0)=0$. With this input, Algorithm
3 returns

$y(\pm\sqrt{3})=0,$

showing that the only non-zero algebraic points where the $E$-function
$J_{0}^{(4)}$ is algebraic are $\pm\sqrt{3}$, where it vanishes. Moreover, the
algorithm described in Section 3.3 provides the canonical decomposition
$J_{0}^{(4)}=p+qg$, where

$p(z)=0,\;q(z)=z^{2}-3\quad\text{and}\quad g(z)=-J_{2}(z)/z^{2}.$ (15)

Here, the purely transcendental $E$-function $g$ (222Siegel first proved that
Bessel’s function $J_{2}$ is purely transcendental in [44, pp. 31-32, §4].) is
given by the differential equation

$zy(z)+5y^{\prime}(z)+y^{\prime\prime}(z)=0,$

with initial conditions $y(0)=-1/8,y^{\prime}(0)=0$. The decomposition (15)
explains the a priori unexpected fact that
$\mathrm{Exc}(J_{0}^{(4)})=\\{\pm\sqrt{3}\\}$.

###  4.3. Two families of $E$-functions

It is easy to construct $E$-functions that take algebraic values at certain
chosen algebraic points: consider $p+qg$ where $p,q\in{\mathbb{Q}}[z]$ and $g$
is any $E$-function in ${\mathbb{Q}}[[z]]$. Conversely, as we have seen in
Theorem 3.4, any $E$-function $f\in{\mathbb{Q}}[[z]]$ can be written $f=p+qg$
where $p,q\in{\mathbb{Q}}[z]$ and $g$ is a purely transcendental $E$-function.

It turns out to be difficult to find an $E$-function which takes an algebraic
value at a non-zero algebraic point and which is not obviously of the form
$p+qg$ as above. The goal of this section is to provide two infinite families
of $E$-functions for which we believe it is difficult to guess a priori
Propositions 4.1 and 4.3 below. Their proof is inspired in part by that of the
evaluation $J_{0}^{(4)}(\pm\sqrt{3})=0$ above. Besides their theoretical
interest, we used these propositions to check the correctness of various
routines of our algorithms.

We start with a result on the exceptional values of an infinite family of
${}_{1}F_{1}$ functions.

######  Proposition 4.1.

Let $a\in\mathbb{Q}\setminus\mathbb{Z}_{\leq 0}$ and
$d\in\mathbb{N}\setminus\\{0\\}$. Then:

  1. $(i)$

$R(z):=\sum_{k=0}^{d}\binom{d}{k}(a)_{k}z^{d-k}$ has $d$ simple roots;

  2. $(ii)$

$\mathrm{Exc}({}_{1}F_{1}[d+1;a+d+1;-z])$ coincides with the set of roots of
$R$;

  3. $(iii)$

for any root $\rho$ of $R$, the following identity holds:

${}_{1}F_{1}[d+1;a+d+1;-\rho]=-\frac{(a)_{d+1}}{\rho R^{\prime}(\rho)}.$ (16)

###### Proof.

From the differential equation
$zy^{\prime\prime}(z)+(a+z+d+1)y^{\prime}(z)+(d+1)y(z)=0$ satisfied by
${}_{1}F_{1}[d+1;a+d+1;-z]$, Algorithm 2 computes its adjoint
$zy^{\prime\prime}(z)-(z+a+d-1)y^{\prime}(z)+dy(z)=0$ and discovers that it
admits $R$ as a non-zero polynomial solution. From there, it computes
$b_{1}=zR(z)$ and $b_{0}=(a+z+d+1)R(z)-b_{1}^{\prime}(z)$, with the property
that ${}_{1}F_{1}[d+1;a+d+1;-z]$ is a solution of the inhomogeneous
differential equation

$b_{1}(z)y^{\prime}(z)+b_{0}(z)y(z)=a(a+1)\cdots(a+d).$ (17)

It follows from (17) that $R$ has only simple roots, since if
$R(\rho)=R^{\prime}(\rho)=0$, then $b_{1}(\rho)=b_{0}(\rho)=0$, hence
$a\in\\{0,-1,\ldots,-d\\}$, which is impossible. This proves $(i)$.

Next, Algorithm 3 evaluates (17) at the roots $\rho$ of $R$. Since
$b_{1}(\rho)=0$, it follows that
${}_{1}F_{1}[d+1;a+d+1;-\rho]=a(a+1)\cdots(a+d)/b_{0}(\rho)$, and hence
$(iii)$ holds.

To prove $(ii)$, note that $f(z):={}_{1}F_{1}[d+1;a+d+1;-z]$ is a
transcendental $E$-function such that $1,f,f^{\prime}$ are linearly dependent
over $\mathbb{Q}(z)$ (as Eq. (17) shows), and $(f,f^{\prime})^{\mathrm{T}}$ is
solution of a differential system with only 0 as singularity. In particular,
by the Siegel-Shidlovskii theorem, for any
$\alpha\in{\overline{\mathbb{Q}}}^{*}$ such that
$f(\alpha)\in{\overline{\mathbb{Q}}}$ we have
$f^{\prime}(\alpha)\notin{\overline{\mathbb{Q}}}$, and consequently the
differential equation (17) shows that
$\mathrm{Exc}({}_{1}F_{1}[d+1;a+d+1;-z])$ coincides with the set of roots of
$R$, proving $(ii)$. ∎

######  Remark 4.2.

_Let us now make a few remarks on Proposition 4.1._

  1. (1)

We do not know if the evaluation (16) is available in the (very rich)
literature on special functions. It is remarkable that it was discovered (and
proved) using our algorithms. Note that Eq. (16) holds more generally for all
$a\in\mathbb{C}\setminus\mathbb{Z}_{\leq 0}$.

  2. (2)

Parts $(i)$ and $(ii)$ of Proposition 4.1 also hold when $d=0$, in the sense
that $\mathrm{Exc}({}_{1}F_{1}[1;a+1;-z])=\emptyset$ for all
$a\in\mathbb{Q}\setminus\mathbb{Z}_{\leq 0}$, see e.g. [43, p. 185, Theorem
1].

  3. (3)

The polynomial $R(z)$ in Proposition 4.1 is equal to
$(a)_{d}\cdot{}_{1}F_{1}[-d;1-a-d;z]$, thus it can be expressed in terms of
generalized Laguerre polynomials as $R(z)=(-1)^{d}d!\cdot L_{d}^{(-a-d)}(z)$.
As proved by Schur [42], the discriminant of $R$ is equal to
$\prod_{j=2}^{d}j^{j}(j-a-d)^{j-1}$, which is non-zero since
$a\notin{\mathbb{Z}}_{\leq 0}$; this yields a different proof that $R$ has
only single roots.

  4. (4)

When $d=1$, the rational canonical decomposition of
$f(z):={}_{1}F_{1}[2;a+2;-z]$ given by Theorem 3.4 is $f=p+qg$ with $p=a+1$,
$q=z+a$ and $g=-{}_{1}F_{1}[1;a+2;-z]$ (note that $g$ is purely transcendental
by remark (2) above).

Decompositions of $f(z):={}_{1}F_{1}[d+1;a+d+1;-z]$ can easily be written down
when $d\geq 2$ but they are neither as explicit nor necessarily canonical.
Since $b_{0}$ and $b_{1}$ (in the proof of Proposition 4.1) are coprime, there
exist $u,v\in\mathbb{Q}[z]$ such that $b_{1}u+b_{0}v=1$. Then we have the
decomposition $f=(a)_{d+1}v+Rg$, where $g(z):=z(u(z)f(z)-v(z)f^{\prime}(z))$
is purely transcendental. Indeed, the decomposition is immediate to check and
let $\alpha\in{\overline{\mathbb{Q}}}^{*}$ be such that
$g(\alpha)\in{\overline{\mathbb{Q}}}$. Then
$f(\alpha)\in{\overline{\mathbb{Q}}}$ as well, hence $b_{1}(\alpha)=0$ by
$(ii)$ in Proposition 4.1, so that $v(\alpha)\neq 0$ by the relation
$b_{1}u+b_{0}v=1$. Therefore $g(\alpha)\notin{\overline{\mathbb{Q}}}$ because
$f^{\prime}(\alpha)\notin{\overline{\mathbb{Q}}}$. This contradiction proves
that there is no such $\alpha$.

When $d=1$, this procedure provides an alternative way to obtain the above
canonical decomposition of $f(z):={}_{1}F_{1}[2;a+2;-z]$, with $g$ represented
as $g(z)=z/(a(a+1))\cdot((z+a-1)f(z)+(z-1)f^{\prime}(z))-1$.

When $d=2$, it provides a decomposition of $f(z):={}_{1}F_{1}[3;a+3;-z]$ as
$f=p+qg$, where $p(z)=z^{2}/2+(a-2)z/2+1$, $q(z)=z^{2}+2az+a(a+1)$ and
$g(z)=-z/(2a(a+1)(a+2))\cdot((z^{2}+2(a-1)z+a^{2}+2)f+(z^{2}+(a-2)z+2)f^{\prime})$.
The canonical decomposition of $f$ is then readily obtained as
$f=\tilde{p}+q\tilde{g}$, where $\tilde{p}:=p-q/2$ and $\tilde{g}:=g+1/2$.

  5. (5)

When $d=2$, Proposition 4.1 implies the following evaluation:

$e^{a\mp i\sqrt{a}}{}_{1}F_{1}[a;a+3;-a\pm i\sqrt{a}]={}_{1}F_{1}[3;a+3;a\mp
i\sqrt{a}]=(a+2)(1\mp i\sqrt{a})/2.$ (18)

The left-hand side is a special case of Kummer’s identity
$e^{-z}{}_{1}F_{1}[a,b,z]={}_{1}F_{1}[b-a,b,-z]$. The right-hand side follows
from the fact that the roots of $R(z)=z^{2}+2az+a(a+1)$ are $\\{-a\pm
i\sqrt{a}\\}$, since then Proposition 4.1 implies
${}_{1}F_{1}[3;a+3;\rho]=-a(a+1)(a+2)/(\rho
R^{\prime}(\rho))=(a+1)(a+2)/(2(\rho+a+1))$.

  6. (6)

Generalized Laguerre polynomials are most of the time irreducible in
${\mathbb{Q}}[z]$, but not always. Filaseta and Lam [21, Thm. 1] proved that
if $\alpha\notin{\mathbb{Z}}_{<0}$, then $L_{d}^{(\alpha)}(z)$ is irreducible
in ${\mathbb{Q}}[z]$ for sufficiently large $d$. However, for some values
$d,\alpha$, the polynomial $L_{d}^{(\alpha)}(z)$ can be reducible. This is so
e.g. for $d=5,\alpha=7/5$, for which $L_{d}^{(\alpha)}(z)$ admits the linear
factor $z-12/5$. This observation leads to simple particular cases of (16)
such as:

$\displaystyle{}_{1}F_{1}[6;-2/5;-12/5]$ $\displaystyle=1309/625,$ (19)
$\displaystyle{}_{1}F_{1}[6;-34/5;-84/5]$ $\displaystyle=8437/625,$ (20)
$\displaystyle{}_{1}F_{1}[5;-113/3;-140/3]$ $\displaystyle=-30073/27,$ (21)
$\displaystyle{}_{1}F_{1}[4;1/7;-6/7]$ $\displaystyle=-65/49,$ (22)
$\displaystyle{}_{1}F_{1}[4;703/725;-312/725]$ $\displaystyle=-20999/525625.$
(23)

Note that such “rational evaluations” are quite rare. We systematically
searched for those arising from Laguerre polynomials $L_{d}^{(\alpha)}$ with
$\alpha\in{\mathbb{Q}}_{\geq 0}$ having numerator and denominator bounded by
$1000$. With $d=3$, we could find 46 such identities, of which (22) is the
simplest and (23) is one of the most complicated. With $d=4$, we could only
find one such identity, namely (21), while with $d=5$ we found the two
identities (19) and (20). We did not find any such identity with $d>5$. For
$d=2$, all rational identities that we have found belong to the following
infinite family:

${}_{1}F_{1}[3;11/4-n^{2}-n;-(2n+3)(2n+1)/4]=(2n-1)(4n^{2}+4n-7)/16.$ (24)

Note that this identity is a particular case of (18), and holds for all
$n\in{\mathbb{Q}}$ for which the left-hand side of (24) is well-defined.

Classifying all pairs $(d,a)\in\mathbb{N}\times{\mathbb{Q}}$ such that the
$E$-function ${}_{1}F_{1}[d+1;a+d+1;-z]$ takes algebraic values at _rational_
points $z$ is a non-trivial task, since by Proposition 4.1 this is equivalent
to finding $(d,a)\in\mathbb{N}\times{\mathbb{Q}}$ such that $L_{d}^{(-a-d)}$
admits a rational root.

  7. (7)

With other choices such as $d=4$, $\alpha=12/5$, the polynomial
$L_{d}^{(\alpha)}$ has quadratic factors. In particular, we obtain the
evaluation

${}_{1}F_{1}[5;-7/5;(6\sqrt{15}-42)/5]={11}/{5}+{66\sqrt{15}}/{125},$ (25)

which is not a particular case of (18).

We were unable to locate in previous works any of the identities (19)–(25),
including in online encyclopedias such as Wolfram’s mathematical functions
site and the Digital Library of Mathematical Functions. Given how vast the
literature on special functions is, we would not be surprized that some of
these identities were already tabulated.

  8. (8)

Equation 16 can also be proved directly starting from the relation between
these ${}_{1}F_{1}$ and the incomplete gamma function [36, 13.6.5], [3,
13.6.10]:

$f(z):=\frac{1}{a}{}_{1}F_{1}[1,a+1,-z]=(-z)^{-a}e^{-z}\gamma(a,-z).$

Successive differentiation of the hypergeometric series shows that

$f^{(d)}(z)=\frac{(-1)^{d}d!}{(a)_{d}}{}_{1}F_{1}[d+1,a+d+1,-z].$

On the other hand, we have
$\left(\gamma(a,-z)\right)^{\prime}=(-z)^{a}e^{z}/z$ by [36, 8.1], [3, 6.5.2].
Thus, by induction, there are two families of polynomials $(R_{d})$ and
$(Q_{d})$ such that

$\frac{(-1)^{d}d!}{(a)_{d}}{}_{1}F_{1}[d+1,a+d+1,-z]=R_{d}(z)e^{-z}(-z)^{-a-d}\gamma(a,-z)+(-z)^{-d}Q_{d}(z)$
(26)

with

$R_{d}(z)=e^{z}(-z)^{a+d}\left((-z)^{-a}e^{-z}\right)^{(d)},\quad
Q_{d+1}=dQ_{d}-zQ_{d}^{\prime}+R_{d},\quad Q_{0}=0.$

The polynomial $R_{d}$ is exactly the Laguerre polynomial $R$ from before.

Thus, Eq. 16 boils down to an evaluation of this more general formula at a
root $\rho$ of $R_{d}$, giving

${}_{1}F_{1}[d+1,a+d+1,-\rho]=\frac{(a)_{d}Q_{d}(\rho)}{d!\rho^{d}}.$

This is not exactly the same formula as above. The proof is completed by
proving by induction that both $Q_{d}$ and $R_{d}$ satisfy the same recurrence
$u_{d+1}=(z+d+a)u_{d}-zdu_{d-1}$, giving an explicit evaluation of the
determinant

$\left|\begin{matrix}Q_{d+1}&R_{d+1}\\\
Q_{d}&R_{d}\end{matrix}\right|=\left|\begin{matrix}z+d+a&-zd\\\
1&0\end{matrix}\right|\cdot\left|\begin{matrix}Q_{d}&R_{d}\\\
Q_{d-1}&R_{d-1}\end{matrix}\right|=\dots=d!z^{d}.$

Evaluating at $z=\rho$ gives $Q_{d}(\rho)/(d!\rho^{d})=-1/R_{d+1}(\rho)$. This
gives another simple expression for the right-hand side of Eq. 16, which
follows from $R_{d+1}(z)=(z+d+a)R_{d}(z)-zR_{d}^{\prime}(z)$.

The next result considers exceptional values of second derivatives of products
of ${}_{1}F_{1}$ with the exponential function. We recall that
$J_{0}(-iz/2)=e^{-z/2}{}_{1}F_{1}[1/2;1;z]$ is such a product.

######  Proposition 4.3.

Let $c\in{\overline{\mathbb{Q}}}^{*}$ and
$a,b\in\mathbb{Q}\setminus\mathbb{Z}_{\leq 0}$ with $a-b\notin\mathbb{N}$. Let
$F(z):=e^{-cz}{}_{1}F_{1}[a;b;z]$. Then
$\mathrm{Exc}(F^{\prime\prime})=\emptyset$, except in the following (disjoint)
cases:

  1. 1.

if $b=a(2c-1)/c^{2}$, then $\mathrm{Exc}(F^{\prime\prime})=\\{-a/c^{2}\\}$ and
$F^{\prime\prime}(-a/c^{2})=0$;

  2. 2.

if $c=1$ and $b=a+1$, then $\mathrm{Exc}(F^{\prime\prime})=\\{-a\pm
i\sqrt{a}\\}$ and $F^{\prime\prime}(-a\pm i\sqrt{a})=1/(1\pm i\sqrt{a})$.

###### Proof.

With the assumptions on $a,b,c$, we prove below that:

Fact 1. $F$ and $F^{\prime}$ are linearly independent over
${\overline{\mathbb{Q}}}(z)$, i.e. $F$ does not satisfy any homogeneous ODE of
order less than $2$.

Fact 2. $F$ satisfies the second-order ODE

$zF^{\prime\prime}(z)=(z-2cz-b)F^{\prime}(z)+(cz+a-c^{2}z-cb)F(z),$ (27)

and it does not satisfy any inhomogeneous ODE of order 1, unless $c=1$ and
$b=a+1$.

Postponing for a moment the proof of these facts, we distinguish two cases.

Case 1. We first assume that either $c\neq 1$ or $b\neq a+1$. By Fact 1, the
function $F$ is a non-polynomial $E$-function (hence a transcendental one). By
Fact 2 and by Beukers’ Corollary 1.4 of [14], it follows that the numbers $1$,
$F(\xi)$ and $F^{\prime}(\xi)$ are linearly independent over
${\overline{\mathbb{Q}}}$ for any $\xi\in{\overline{\mathbb{Q}}}^{*}$.

Assume now that $b\neq(2ac-a)/c^{2}$. Then, for any $\xi\in\mathbb{C}$, the
numbers $\xi-2c\xi-b$ and $c\xi+a-\xi c^{2}-cb$ cannot be simultaneously equal
to 0. Since the numbers $1$, $F(\xi)$ and $F^{\prime}(\xi)$ are linearly
independent over ${\overline{\mathbb{Q}}}$ for any
$\xi\in{\overline{\mathbb{Q}}}^{*}$, it follows that

$F^{\prime\prime}(\xi)=\frac{1}{\xi}\big{(}(\xi-2c\xi-b)F^{\prime}(\xi)+(c\xi+a-\xi
c^{2}-cb)F(\xi)\big{)}\notin{\overline{\mathbb{Q}}}.$

Hence $\mathrm{Exc}(F^{\prime\prime})=\emptyset$ when $b\neq(2ac-a)/c^{2}$.

It remains to treat the sub-case $b:=(2ac-a)/c^{2}$. With $z:=-a/c^{2}$, we
see that

$z-2cz-b=cz+a-c^{2}z-cb=0,$

so that $F^{\prime\prime}(-a/c^{2})=0$. Since $z-2cz-b$ and $cz+a-c^{2}z-cb$
vanish simultaneously for no other value of $z$ and since the numbers $1$,
$F(\xi)$ and $F^{\prime}(\xi)$ are linearly independent over
${\overline{\mathbb{Q}}}$ for any
$\xi\in{\overline{\mathbb{Q}}}\setminus\\{0,-a/c^{2}\\}$, we deduce that
$F^{\prime\prime}(\xi)\notin{\overline{\mathbb{Q}}}$ for such $\xi$. Hence
$\mathrm{Exc}(F^{\prime\prime})=\\{-a/c^{2}\\}$ when $b=(2ac-a)/c^{2}$. The
proposition is thus proved in Case 1.

Case 2. We assume that $c=1$ and $b=a+1$. Then,
$F^{\prime\prime}(z)=\frac{2}{(a+1)(a+2)}e^{-z}{}_{1}F_{1}[a;a+3;z]$, so
$\mathrm{Exc}(F^{\prime\prime})=\mathrm{Exc}(e^{-z}{}_{1}F_{1}[a;a+3;z])=\\{-a\pm
i\sqrt{a}\\}$. The last equality is a consequence of Proposition 4.1 with
$d=2$. The equality $F^{\prime\prime}(-a\pm i\sqrt{a})=1/(1\pm i\sqrt{a})$
follows from (18).

Proof of Fact 1. When $a,b\in\mathbb{Q}\setminus\mathbb{Z}_{\leq 0}$ and
$a-b\notin\mathbb{N}$, the asymptotic behavior of ${}_{1}F_{1}[a;b;z]$ as
$z\to\infty$ (with $-\pi/2<\arg(z)<3\pi/2$) is given in [3, p. 508, Eq.
13.5.1], [36, 13.7.2]. In the particular cases $z\to\pm\infty$, it reads

${}_{1}F_{1}[a;b;z]\sim_{z\to+\infty}\frac{\Gamma(b)}{\Gamma(a)}e^{z}z^{a-b}\quad\text{and}\quad{}_{1}F_{1}[a;b;z]\sim_{z\to-\infty}\frac{\Gamma(b)}{\Gamma(b-a)}e^{i\pi
a}z^{a}.$

This rules out the possibility that $e^{-cz}{}_{1}F_{1}[a;b;z]$ satisfies a
differential equation of order 1 over ${\overline{\mathbb{Q}}}(z)$.

Note that a different (purely algebraic) proof is possible, based on a
reasoning similar to the one in the statement and proof of [15, Lemma 4.2].

Proof of Fact 2. Since ${}_{1}F_{1}[a;b;z]$ satisfies
$zy^{\prime\prime}(z)+(b-z)y^{\prime}(z)-ay(z)=0$ and $e^{-cz}$ satisfies
$y^{\prime}(z)+cy(z)=0$, it follows by a simple computation that $F$ satisfies
(27). By Fact 1, (27) is the minimal-order homogeneous ODE satisfied by $F$.

Assume now that $F$ satisfies an inhomogeneous ODE of order $1$. We will
follow the reasoning in Algorithm 2, and show that the adjoint of (27) does
not possess any non-zero rational solutions in ${\overline{\mathbb{Q}}}(z)$,
unless $c=1$ and $b=a+1$.

The adjoint equation of (27) writes

$zy^{\prime\prime}(z)+((1-2c)z-b+2)y^{\prime}(z)+(c(c-1)z+bc-a-2c+1)y(z)=0.$
(28)

If it admits a rational solution $R(z)\in{\overline{\mathbb{Q}}}(z)$, then the
only potential pole of $R$ can be located at $z=0$. The indicial equation of
(28) at $z=0$ is $s(s-b+1)$. Hence, the possible valuations at $z=0$ of $R$
are $0$ and $b-1$. Since $b\in\mathbb{Q}\setminus\mathbb{Z}_{\leq 0}$, this
implies that $R$ is actually a polynomial solution in
${\overline{\mathbb{Q}}}[z]$ of (28). If $c\neq 1$, then the indicial
polynomial at infinity of (28) is a non-zero constant, equal to $c^{2}-c$;
therefore in that case, $R$ cannot be a polynomial solution. It follows that
$c=1$. Now, the indicial polynomial at infinity of (28) is $s-a+b-1$, hence
the only possible degree of $R$ is $a-b+1$. Since $a-b\notin\mathbb{N}$, this
implies that $b=a+1$ and that $R$ is a constant in ${\overline{\mathbb{Q}}}$.
In this case, (28) admits the rational solution $y(z)=1$, and $F$ satisfies
$zF^{\prime}(z)+(z+a)F(z)=a$. ∎

Note that, in the spirit of Proposition 4.3, the following examples that can
be treated along the same lines: the third derivatives of

$e^{-z/9}{}_{2}F_{2}[1/144,1/144;-7/16,-7/16;z]\quad\textup{and}\quad
e^{-z/3}{}_{2}F_{2}[1/4,3/4;5/4,-9/4;z]$

vanish at $z=3/16$ and $z=-9/4$, respectively. Our minimization algorithm
finds their minimal differential equations, which are too big to be written
here: their (order, degree) are $(3,8)$ and $(2,7)$, respectively. Algorithm 3
then shows that they are transcendental functions and that $\\{3/16\\}$ and
$\\{-9/4\\}$ are the exceptional values sets in each case. A complete
classification as in Proposition 4.3 seems to be currently out of reach
though. Identities like ${}_{2}F_{2}[1/2,1/3;-1/2,-2/3;z]=(1-z/2+3z^{2})e^{z}$
show that an assumption corresponding to the assumption $a-b\notin\mathbb{N}$
in Proposition 4.3 is obviously necessary on the rational parameters of the
${}_{2}F_{2}[a,b;c,d;z]$ to avoid trivial situations: besides the fact that
$a,b,c,d$ must all be in $\mathbb{Q}\setminus\mathbb{Z}_{\leq 0}$, we must not
have $a-c\in\mathbb{N}$ and $b-d\in\mathbb{N}$, or $a-d\in\mathbb{N}$ and
$b-c\in\mathbb{N}$ (note that in the second example above, we have
$3/4-(-9/4)\in\mathbb{N}$ but $1/4-5/4\notin\mathbb{N}$, while
$1/4-(-9/4)\notin\mathbb{N}$ and $3/4-5/4\notin\mathbb{N}$).

###  4.4. An example with Gauss’ hypergeometric function

The approach leading to special evaluations is very general and not restricted
to $E$-functions. For instance, Gauss’ hypergeometric function
${}_{2}F_{1}[a,b;c;z]$ satisfies a differential equation whose adjoint is
solved by $R(z):={}_{2}F_{1}[1-a,1-b;2-c;z]$. The approach from Section 2.4
then deduces that the hypergeometric function satisfies a first order
_inhomogeneous_ equation, with coefficients that are not polynomials in
general, namely

$z(z-1)R(z)y^{\prime}(z)+(z(1-z)R^{\prime}(z)+((a+b-1)z+1-c)R(z))y(z)+c-1=0.$

It follows that if $\rho$ is a simple zero of $R(z)$ different from 0,1, one
gets the special evaluation

${}_{2}F_{1}[a,b;c;\rho]=\frac{1-c}{\rho(1-\rho)R^{\prime}(\rho)}.$

The special case $c=a+k+1$ ($k\in\mathbb{N}$) gives a nice analogue of
Proposition 4.1. To state it, recall that the $k$-th _Jacobi polynomial_
$P_{k}^{(\alpha,\beta)}$ with parameters $\alpha,\beta\in\mathbb{C}$ is
defined by

$P_{k}^{(\alpha,\beta)}(z):=2^{-k}\cdot\sum_{j=0}^{k}\binom{k+\alpha}{k-j}\binom{k+\beta}{j}(z-1)^{j}(z+1)^{k-j}.$

It is classical [46, §6.72] that $P_{k}^{(\alpha,\beta)}$ has only simple
roots (which are even real and in the interval $(-1,1)$ if $\alpha$ and
$\beta$ are both real and greater than $-1$), with the notable exception of
$\pm 1$ which is a multiple root of $P_{k}^{(\alpha,\beta)}$ if one of the
parameters $\alpha,\beta$ is in $\\{-1,\ldots,-k\\}$.

######  Proposition 4.4.

Let $a\in\mathbb{Q}\setminus\mathbb{Z}_{\leq 0}$, $b\in\mathbb{Q}$ and
$k\in\mathbb{N}$. If $\rho\in{\overline{\mathbb{Q}}}\setminus\\{0,1\\}$ is a
root of the polynomial $P_{k}^{(-k-a,b-k-1)}(1-2z)$, then

${}_{2}F_{1}[a,b;a+k+1;\rho]=\frac{(-1)^{k}a\binom{a+k}{k}(1-\rho)^{k-b}}{(k+a-b)P_{k}^{(-k-a,b-k)}(1-2\rho)}.$
(29)

###### Proof.

With $c=a+k+1$, the hypergeometric function $R(z)$ becomes

$R(z)={}_{2}F_{1}[1-a,1-b;1-a-k;z]=\frac{(-1)^{k}k!}{(a)_{k}}{(1-z)^{b-k-1}P_{k}^{(-k-a,b-k-1)}(1-2z)}.$

Its roots different from 0,1 are the roots of $P_{k}^{(-k-a,b-k-1)}(z)$
different from $-1,1$, which are all simple. The formula for the denominator
comes from the derivative

$R^{\prime}(z)=-\frac{(-1)^{k}}{z\binom{a+k-1}{k}}\left((k+a-b)(1-z)^{b-k-1}P_{k}^{(-k-a,b-k)}(1-2z)\right.\\\
\left.+((a-1)z+b-a-k)(1-z)^{b-k-2}P_{k}^{(-k-a,b-k-1)}(1-2z)\right).\qed$

######  Remark 4.5.

_Let us conclude with a few remarks on Proposition 4.4._

  1. (1)

If none of $a,b,a-b$ is an integer, then $f(z)={}_{2}F_{1}[a,b;a+k+1;z]$ is a
transcendental function [53]. Therefore, the evaluation in Eq. 29 provides
very simple particular cases of algebraic values taken by transcendental
$G$-functions at algebraic points. Note that (29) holds more generally for
$a\in\mathbb{C}\setminus\mathbb{Z}_{\leq 0}$ and $b\in\mathbb{C}$.

  2. (2)

As in the case of Proposition 4.1, another proof of Proposition 4.4 relies on
a relation analogous to Eq. 26 between ${}_{2}F_{1}[a,b;a+k+1;z]$,
${}_{2}F_{1}[a,b;a;z]=(1-z)^{-b}$ and

${}_{2}F_{1}[a,b;a+1;z]=az^{-a}B_{z}(a,1-b)=az^{-a}\int_{0}^{z}{z^{a-1}(1-z)^{b}\,dz},$

an incomplete beta function. The relation is obtained from those two by
repeated use of a contiguity relation.

  3. (3)

As in Remark 4.2, nice special cases of (29) can be obtained by studying
triples $(a,b,k)$ for which the Jacobi polynomial $P_{k}^{(-k-a,b-k-1)}$
factors non-trivially. For instance, the triples $(2/5,3/5,5)$ and
$(2/3,7/3,3)$ yield the evaluations

$\displaystyle{}_{2}F_{1}\left[\frac{2}{5},\frac{3}{5};\frac{32}{5};\frac{1}{2}\right]$
$\displaystyle=\frac{1683}{2500}\,\sqrt[5]{8}$ and
$\displaystyle{}_{2}F_{1}\left[\frac{2}{3},\frac{7}{3};\frac{14}{3};\frac{3\sqrt{5}-5}{2}\right]$
$\displaystyle=\frac{44}{27\sqrt[3]{28-12\sqrt{5}}}.$

(Here, by item (1) above, both hypergeometric functions are transcendental.)

###  4.5. Implementation issues

The main issue is to avoid the computation of high-order expansions of power
series with rational coefficients. A first gain is achieved by performing most
of the computation modulo a sufficiently large prime number (we take a 31-bit
long prime number). When a factor is found with modular coefficients, then the
actual degree bounds from that factor are used to determine how many rational
coefficients of the power series have to be computed and then obtain the
differential operator with rational coefficients. Another gain is achieved by
not necessarily computing expansions at orders up to the bound obtained from
the degree bounds. Instead, one keeps doubling the number of coefficients
computed previously until the bound is reached or until a factor is found.
This induces important savings in the latter case and not too large an
additional cost otherwise.

Another place where some time can be saved is in the optimization problems.
The computation of an approximant basis returns a linear differential operator
of small order if one exists with the given degree bounds. Thus the
computation of a tight bound on the number of apparent singularities is only
useful if it leads to a bound on the degrees smaller than a previously known
one. One can therefore add an extra inequality to the optimization problem so
that the solver does not waste time computing an optimum which is larger than
what is already known.

###  4.6. Timings

$(m,p)$ | (ord,deg) | (ord,deg) | (ord,deg) | number | number | time | time
---|---|---|---|---|---|---|---
| rec | original | minimal | modular | rational | minim. | alg. values
|  | diff.eq. | diff.eq. | terms | terms | (s.) | (s.)
(1,1) | (2,2) | (2,1) | (2,1) | – | – | 0.12 | 0.002
(1,2) | (3,5) | (5,3) | (4,3) | 47 | 30 | 0.21 | 0.01
(1,3) | (4,10) | (10,8) | (6,8) | 97 | 75 | 0.46 | 0.04
(1,4) | (5,19) | (19,16) | (9,16) | 706 | 185 | 3.4 | 0.2
(2,1) | (2,3) | (3,2) | (3,2) | – | – | 0.04 | 0.001
(2,2) | (2,4) | (4,3) | (4,3) | – | – | 0.05 | 0.002
(2,3) | (4,14) | (14,12) | (8,12) | 326 | 131 | 1.7 | 0.08
(2,4) | (5,26) | (26,24) | (11,24) | 1501 | 317 | 11. | 0.8
(3,1) | (4,10) | (10,8) | (6,8) | 97 | 75 | 0.44 | 0.04
(3,2) | (5,19) | (19,16) | (9,16) | 706 | 185 | 3.1 | 0.2
(3,3) | (6,28) | (28,25) | (12,29) | 2461 | 408 | 22. | 1.4
(3,4) | (7,51) | (51,47) | (16,47) | 1261 | 838 | 123. | 46.
(4,1) | (4,16) | (16,13) | (9,13) | 325 | 155 | 2.9 | 0.2
(4,2) | (5,27) | (27,24) | (12,24) | 1927 | 343 | 16. | 0.6
(4,3) | (6,41) | (41,38) | (15,38) | 5362 | 645 | 73. | 9.1
(4,4) | (6,46) | (46,43) | (18,47) | 9634 | 936 | 212. | 26.
Table 1. Experimental results

Experimental results (333The timings were obtained with Maple2021 on a 2018
Mac mini.) on the family of power series

$f_{m,p}(z)=\sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}\binom{n}{k}^{m}\binom{n+k}{k}^{p}\right)\frac{z^{n}}{n!}$

are reported in Table 1. These power series are exponential generating
functions of Apéry-like sequences, hence they are $E$-functions by design. The
case $(m,p)=(2,1)$ was considered by Adamczewski and Rivoal in [5, p. 706],
who proved that $f_{2,1}$ is a purely transcendental $E$-function. We used our
algorithms to reprove this result and to extend it to other values of $m$ and
$p$, see Table 1.

For each $(m,p)$, we indicate the order and the degree of the (minimal)
recurrence computed by Zeilberger’s algorithm, the order and the degree of the
differential equation deduced from this recurrence, and those of the minimal-
order (homogeneous) differential equation obtained by our implementation when
run on this differential equation. We also give the number of coefficients of
the sequence that were computed modulo a prime number and the number of
rational coefficients of the sequence that were computed. These numbers are
replaced by ‘–’ when no term was computed, minimality having been concluded
directly from consideration of degree bounds. The next column contains the
time in seconds spent in the mininization process (including searching for a
minimal inhomogeneous equation). Finally, the last column gives the extra time
required to search for the exceptional algebraic values for these functions
(and prove that there are none).

In more detail, during the computation of the case $(m,p)=(4,4)$, the steps
that take more that 3 seconds are: 132 sec. for the computation of modular
approximant bases; 54 sec. for the computation of a rational approximant
basis; 9 sec. for the computation of rational terms of the sequence; 6 sec.
for checking that there are no inhomogeneous equations of smaller order; 26
sec. in Beukers’ algorithm (the matrix has dimension 19 and the apparent
singularity has degree 32).

On the basis of these experimental results, we ask the following questions on
the family $f_{m,p}(z)$ and leave them for further research. The data in Table
1, plus a few more experiments (not included in Table 1), are in favor of
positive answers to all these questions.

######  Question 4.6.

Is the $E$-function $f_{m,p}$ purely transcendental for any $m\geq 1$ and
$p\geq 1$?

######  Question 4.7.

Is it true that for any odd $m$, the minimal-order linear differential
equation $L^{\textrm{min}}_{m,p}(y)=0$ satisfied by $f_{m,p}$ has order
$\lfloor(N+1)^{2}/4\rfloor$ and degree $\lfloor N(2N^{2}-3N+4)/12\rfloor$,
where $N=m+p$? In particular, when both $m$ and $p$ are odd, is it true that
${\operatorname{ord}}(L^{\textrm{min}}_{m,p})={\operatorname{ord}}(L^{\textrm{min}}_{p,m})$
and $\deg(L^{\textrm{min}}_{m,p})=\deg(L^{\textrm{min}}_{p,m})$?

## Acknowledgements.

This work has been supported by the French ANR grant DeRerumNatura,
ANR-19-CE40-0018. We thank Frits Beukers and Lucia Di Vizio for many
interesting discussions.

## References

  * [1] Abramov, S. A. Rational solutions of linear differential and difference equations with polynomial coefficients.  Zh. Vychisl. Mat. i Mat. Fiz. 29, 11 (1989), 1611–1620, 1757.
  * [2] Abramov, S. A., and Kvansenko, K. Y. Fast algorithms for the search of the rational solutions of linear differential equations with polynomial coefficients.  In ISSAC’91—Proceedings of the 16th International Symposium on Symbolic and Algebraic Computation (1991), ACM Press, pp. 267–270.
  * [3] Abramowitz, M., and Stegun, I. A. Handbook of mathematical functions with formulas, graphs, and mathematical tables.  National Bureau of Standards Applied Mathematics Series, No. 55. U. S. Government Printing Office, Washington, D.C., 1964.
  * [4] Adamczewski, B., and Faverjon, C. Mahler method: linear relations, transcendence and applications to automatic numbers.  Proc. Lond. Math. Soc. (3) 115, 1 (2017), 55–90.
  * [5] Adamczewski, B., and Rivoal, T. Exceptional values of $E$-functions at algebraic points.  Bull. Lond. Math. Soc. 50, 4 (2018), 697–708.
  * [6] André, Y. Séries Gevrey de type arithmétique. I. Théorèmes de pureté et de dualité.  Ann. of Math. (2) 151, 2 (2000), 705–740.
  * [7] André, Y. Séries Gevrey de type arithmétique. II. Transcendance sans transcendance.  Ann. of Math. (2) 151, 2 (2000), 741–756.
  * [8] André, Y. Solution algebras of differential equations and quasi-homogeneous varieties: a new differential Galois correspondence.  Ann. Sci. Éc. Norm. Supér. (4) 47, 2 (2014), 449–467.
  * [9] Barkatou, A. A rational version of Moser’s algorithm.  In ISSAC’95—Proceedings of the 20th International Symposium on Symbolic and Algebraic Computation (1995), ACM Press, p. 297–302.
  * [10] Barkatou, M. A., and Maddah, S. S. Removing apparent singularities of systems of linear differential equations with rational function coefficients.  In ISSAC’15—Proceedings of the 40th International Symposium on Symbolic and Algebraic Computation (2015), ACM Press, pp. 53–60.
  * [11] Beckermann, B., Labahn, G., and Villard, G. Shifted normal forms of polynomial matrices.  In ISSAC’99—Proceedings of the 24th International Symposium on Symbolic and Algebraic Computation (1999), ACM Press, pp. 189–196.
  * [12] Bertrand, D. On André’s proof of the Siegel-Shidlovsky theorem.  In Colloque Franco-Japonais: Théorie des Nombres Transcendants (Tokyo, 1998), vol. 27 of Sem. Math. Sci. Keio Univ., Yokohama, 1999, pp. 51–63.
  * [13] Bertrand, D., and Beukers, F. Équations différentielles linéaires et majorations de multiplicités.  Ann. Sci. Éc. Norm. Supér. (4) 18, 1 (1985), 181–192.
  * [14] Beukers, F. A refined version of the Siegel-Shidlovskii theorem.  Ann. of Math. (2) 163, 1 (2006), 369–379.
  * [15] Beukers, F., Brownawell, W. D., and Heckman, G. Siegel normality.  Ann. of Math. (2) 127, 2 (1988), 279–308.
  * [16] Bostan, A., Lairez, P., and Salvy, B. Multiple binomial sums.  J. Symbolic Comput. 80, part 2 (2017), 351–386.
  * [17] Bostan, A., Rivoal, T., and Salvy, B. Explicit degree bounds for right factors of linear differential operators.  Bull. Lond. Math. Soc. 53, 1 (2021), 53–62.
  * [18] Churchill, R. C., and Kovacic, J. J. Cyclic vectors.  In Differential algebra and related topics (Newark, NJ, 2000). World Sci. Publ., River Edge, NJ, 2002, pp. 191–218.
  * [19] Cope, F. T. Formal Solutions of Irregular Linear Differential Equations. Part II.  Amer. J. Math. 58, 1 (1936), 130–140.
  * [20] Dehghan, A., Assari, S. M., and Shah, M. Gmmcp tracker: Globally optimal generalized maximum multi clique problem for multiple object tracking.  In 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2015), pp. 4091–4099.
  * [21] Filaseta, M., and Lam, T.-Y. On the irreducibility of the generalized Laguerre polynomials.  Acta Arith. 105, 2 (2002), 177–182.
  * [22] Fischler, S., and Rivoal, T. On the values of $G$-functions.  Comment. Math. Helv. 89, 2 (2014), 313–341.
  * [23] Fischler, S., and Rivoal, T. Arithmetic theory of $E$-operators.  J. Éc. Polytech., Math. 3 (2016), 31–65.
  * [24] Fresán, J., and Jossen, P. A non-hypergeometric $E$-function.  Ann. of Math. (2) 194, 3 (2021), 903–942.
  * [25] Gorelov, V. A. Algebraic independence of the values of $E$-functions at singular points and Siegel’s conjecture.  Math. Notes 67, 2 (2000), 138–151.
  * [26] Gorelov, V. A. On the Siegel conjecture for second-order homogeneous linear differential equations.  Math. Notes 75, 4 (2004), 513–529.
  * [27] Grigoriev, D. Y. Complexity of factoring and calculating the GCD of linear ordinary differential operators.  J. Symbolic Comput. 10, 1 (1990), 7–37.
  * [28] Ince, E. L. Ordinary differential equations.  Dover Publications, New York, 1956.  Reprint of the 1926 edition.
  * [29] Jeannerod, C.-P., Neiger, V., and Villard, G. Fast computation of approximant bases in canonical form.  J. Symbolic Comput. 98 (2020), 192–224.
  * [30] Lepetit, G. Le théorème d’André-Chudnovsky-Katz au sens large.  North-W. Eur. J. of Math. 7 (2021), 83–149.
  * [31] Liouville, J. Second mémoire sur la détermination des intégrales dont la valeur est algébrique.  J. Éc. Polytech. 14 (1833), 149–193.
  * [32] Lorch, L., and Muldoon, M. E. Transcendentality of zeros of higher derivatives of functions involving Bessel functions.  Internat. J. Math. Math. Sci. 18, 3 (1995), 551–560.
  * [33] Malgrange, B. Sur la réduction formelle des équations différentielles à singularités irrégulières.  Preprint, 1979.
  * [34] Moser, J. The order of a singularity in Fuchs’ theory.  Math. Z 72 (1959/1960), 379–398.
  * [35] Nesterenko, Y. V., and Shidlovskiĭ, A. B. On the linear independence of values of $E$-functions.  Mat. Sb. 187, 8 (1996), 93–108.
  * [36] Olver, F. W. J., Olde Daalhuis, A. B., Lozier, D. W., Schneider, B. I., Boisvert, R. F., Clark, C. W., Miller, B. R., Saunders, B. V., Cohl, H. S., and M. A. McClain, e. NIST digital library of mathematical functions.  http://dlmf.nist.gov/, 2022.
  * [37] Ore, O. Theory of non-commutative polynomials.  Ann. of Math. 34 (1933), 480–508.
  * [38] Philippon, P. Galois groups and automatic numbers.  J. Lond. Math. Soc., II. Ser. 92, 3 (2015), 596–614.
  * [39] Poole, E. G. C. Introduction to the theory of linear differential equations.  Dover Publications Inc., New York, 1960.
  * [40] Popov, V. M. Invariant description of linear, time-invariant controllable systems.  SIAM J. Control 10 (1972), 252–264.
  * [41] Rivoal, T. Valeurs algébriques de $E$-fonctions aux points algébriques.  Unpublished note (2016), 4 pages.
  * [42] Schur, J. Affektlose Gleichungen in der Theorie der Laguerreschen und Hermiteschen Polynome.  J. Reine Angew. Math. 165 (1931), 52–58.
  * [43] Shidlovskii, A. B. Transcendental numbers, vol. 12 of De Gruyter Studies in Mathematics.  Walter de Gruyter & Co., Berlin, 1989.  Translated from the Russian by Neal Koblitz, with a foreword by W. Dale Brownawell.
  * [44] Siegel, C. L. Über einige Anwendungen diophantischer Approximationen [reprint of Abhandlungen der Preußischen Akademie der Wissenschaften. Physikalisch-mathematische Klasse 1929, Nr. 1].  In On some applications of Diophantine approximations, vol. 2 of Quad./Monogr. Ed. Norm., Pisa, 2014, pp. 81–138.
  * [45] Singer, M. F. Introduction to the Galois theory of linear differential equations.  In Algebraic theory of differential equations, vol. 357 of London Math. Soc. Lecture Note Ser. Cambridge Univ. Press, Cambridge, 2009, pp. 1–82.
  * [46] Szegő, G. Orthogonal polynomials, fourth ed.  American Mathematical Society Colloquium Publications, Vol. XXIII. American Mathematical Society, Providence, R.I., 1975.
  * [47] Van Barel, M., and Bultheel, A. A general module-theoretic framework for vector M-Padé and matrix rational interpolation.  Numer. Algorithms 3, 1-4 (1992), 451–461.
  * [48] van der Hoeven, J. On the complexity of skew arithmetic.  Appl. Algebra Engrg. Comm. Comput. 27, 2 (2016), 105–122.
  * [49] van der Put, M., and Singer, M. F. Galois theory of linear differential equations, vol. 328 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, Berlin, 2003.
  * [50] van Hoeij, M. Factorization of differential operators with rational functions coefficients.  J. Symbolic Comput. 24, 5 (1997), 537–561.
  * [51] van Hoeij, M. Formal solutions and factorization of differential operators with power series coefficients.  J. Symbolic Comput. 24, 1 (1997), 1–30.
  * [52] Vazirani, V. V. Approximation algorithms.  Springer-Verlag, Berlin, 2001.
  * [53] Vidūnas, R. Degenerate Gauss hypergeometric functions.  Kyushu J. Math. 61, 1 (2007), 109–135.
  * [54] Zeilberger, D. The method of creative telescoping.  J. Symbolic Comput. 11, 3 (1991), 195–204.