# First-order factors of linear Mahler operators

Frédéric Chyzak, Thomas Dreyfus, Philippe Dumas, Marc Mezzarobba

► To cite this version:

Frédéric Chyzak, Thomas Dreyfus, Philippe Dumas, Marc Mezzarobba. First-order factors of linear Mahler operators. Journal of Symbolic Computation, 2025, 10.1016/j.jsc.2025.102424. hal-04502052v3

HAL Id: hal-04502052

https://inria.hal.science/hal-04502052v3

Submitted on 3 Feb 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# First-order factors of linear Mahler operators*,**

Frédéric Chyzakª, Thomas Dreyfusᵇ, Philippe Dumasª, Marc Mezzarobbaᶜ

ªInria, France

ᵇUniversité Bourgogne Europe, CNRS, IMB UMR 5584, F-21000, Dijon, France

ᶜLIX, CNRS, École polytechnique, Institut polytechnique de Paris, 91120, Palaiseau, France

## Abstract

We develop and compare two algorithms for computing first-order right-hand factors in the ring of linear Mahler operators $l_r M^r + \dots + l_1 M + l_0$ where $l_0, \dots, l_r$ are polynomials in $x$ and $Mx = x^b M$ for some integer $b \ge 2$. In other words, we give algorithms for finding all formal infinite product solutions of linear functional equations $l_r(x)f(x^{b^r}) + \dots + l_1(x)f(x^b) + l_0(x)f(x) = 0$.

The first of our algorithms is adapted from Petkovšek's classical algorithm for the analogous problem in the case of linear recurrences. The second one proceeds by computing a basis of generalized power series solutions of the functional equation and by using Hermite-Padé approximants to detect those linear combinations of the solutions that correspond to first-order factors.

We present implementations of both algorithms and discuss their use in combination with criteria from the literature to prove the differential transcendence of power series solutions of Mahler equations.

**Keywords:** Mahler operator, factorization, hypergeometric solution, infinite product, Petkovšek's algorithm, Hermite-Padé approximant

## Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>1.1</td><td>Mahler equations</td><td>3</td></tr><tr><td>1.2</td><td>First-order factors and hypergeometric solutions</td><td>4</td></tr></table>

<footnote>*Dedicated to the memory of Marko Petkovšek.</footnote>

<footnote>**CC BY This work is licensed under a Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/).</footnote>

<footnote>For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.</footnote>

<footnote>Supported in part by the French ANR grant De rerum natura (ANR-19-CE40-0018) and by the French-Austrian ANR-FWF grant EAGLES (ANR-22-CE91-0007 & FWF-I-6130-N). This work has been achieved in the frame of the EIPHI Graduate school (contract ANR-17-EURE-0002).</footnote>

<footnote>Email addresses: frederic.chyzak@inria.fr (Frédéric Chyzak), thomas.dreyfus@cnrs.fr (Thomas Dreyfus), philippe.dumas@inria.fr (Philippe Dumas), marc@mezzarobba.net (Marc Mezzarobba)</footnote>
<table><tr><td>1.3 Contribution</td><td>6</td></tr><tr><td>1.4 Related work</td><td>7</td></tr><tr><td>1.5 Outline</td><td>8</td></tr><tr><td>1.6 Notation</td><td>9</td></tr></table>

# I Structural results 9

## 2 Structure of the solution set of the Riccati equation 9

2.1 Basic notions of difference algebra 9

2.2 First-order factors and their solutions 10

2.3 Parametrization of the solution set 12

## 3 Generalized series solutions of the linear equation 14

3.1 The difference ring $\mathcal{D}$ 14

3.2 Bounds 16

3.3 Puiseux-hypergeometric solutions 18

## 4 Ramified rational solutions to the Riccati equation 21

4.1 Hypergeometric elements 21

4.2 Similarity classes 22

## 5 Degree bounds for rational solutions 24

# II Algorithm by exploration of the possible singularities 26

## 6 Mahlerian variant of Petkovšek's algorithm 26

6.1 Petkovšek's classical algorithm 27

6.2 Roques's algorithm for order 2 28

6.3 A new algorithm for higher-order Mahler equations 29

6.4 Existence and computation of bounded Gosper-Petkovšek forms 34

## 7 A more practical algorithm 36

7.1 Ensuring coprimality 36

7.2 Avoiding redundant pairs 38

7.3 Avoiding redundant computations of Newton polygons 41

# III Algorithm by reconstruction from series 43

## 8 Solutions as syzygies 43

8.1 Approximate syzygies 46

8.2 Structured syzygies 48

8.3 A linear-algebra condition on structured syzygies 49

8.4 Candidate solutions 51

8.5 Validating candidate solutions 55

8.6 Ramified rational solutions 56

8.7 A supplementary remark on the rank of syzygy module 57

## 9 Algorithms by Hermite-Padé approximation 58

9.1 Algorithm by sieving candidates for solving the Riccati equation 58

9.2 Rational solving of the linear Mahler equation 64

# IV Implementation and application 66
<table>
  <tbody>
    <tr>
      <td><b>10 Benchmark</b></td>
      <td><b>66</b></td>
    </tr>
    <tr>
      <td>10.1 Implementation</td>
      <td>66</td>
    </tr>
    <tr>
      <td>10.2 Examples</td>
      <td>67</td>
    </tr>
    <tr>
      <td>10.3 Discussion of the timings</td>
      <td>70</td>
    </tr>
    <tr>
      <td><b>11 Differential transcendence of Mahler functions</b></td>
      <td><b>74</b></td>
    </tr>
  </tbody>
</table>

# 1. Introduction

## 1.1. Mahler equations

Mahler equations are a type of functional equations that, for some fixed integer $b \ge 2$, relate the iterates $y(x), y(x^b), y(x^{b^2}), \dots$, of an unknown function $y$ under the substitution of $x^b$ for $x$. Mahler originally considered nonlinear multivariate equations of this type in his work in transcendence theory in the 1920s. The focus later shifted to linear univariate equations, in relation to the study of automatic sequences. More recently, questions in difference Galois theory related to the existence of nonlinear differential equations satisfied by solutions of difference equations of various types led to a revival of the topic.

**Example 1.1.** 1. Let $t_0t_1\cdots = 01101001\cdots$ be the Thue–Morse sequence (OEIS, A010060), defined as the fixed point with $t_0 = 0$ of the substitution $\{0 \mapsto 01, 1 \mapsto 10\}$. It is classical that the series $y(x) = \sum_{n=0}^{\infty} (-1)^{t_n} x^n$ can be written as

$$y(x) = \prod_{j=0}^{\infty} (1 - x^{2^j}).$$

This infinite product obviously satisfies the linear first-order Mahler equation

$$y(x^2) = (1-x)y(x).$$

2. The Baum–Sweet sequence $(a_n)_{n \in \mathbb{N}}$ is the automatic sequence defined by $a_n = 1$ if the binary representation of $n$ contains no blocks of consecutive 0 of odd length, and $a_n = 0$ otherwise (OEIS, A086747). The generating function $y(x) = \sum_{n \in \mathbb{N}} a_n x^n$ satisfies the Mahler equation

$$y(x^4) + xy(x^2) - y(x) = 0.$$

3. The Stern–Brocot sequence $(a_n)_{n \in \mathbb{N}}$ (OEIS, A002487) was introduced by Stern to enumerate the nonnegative rational numbers bijectively by the numbers $a_n/a_{n+1}$. Allouche and Shallit (1992) showed that the sequence is completely defined by $a_0 = 0$, $a_1 = 1$, and for all $n \in \mathbb{N}$,

$$a_{2n} = a_n, \quad a_{4n+1} = a_n + a_{2n+1}, \quad a_{4n+3} = 2a_{2n+1} - a_n.$$

Starting from this property, the algorithms developed in the present article can be used to rediscover the well-known expression

$$y(x) = \sum_{n \ge 0} a_n x^n = x \prod_{k \ge 0} (1 + x^{2^k} + x^{2^{k+1}})$$
as a consequence of the equation

$$xy(x) - (1 + x + 2x^2)y(x^2) + (1 + x^2 + x^4)y(x^4) = 0.$$

4. Adamczewski and Faverjon (2017, Example 8.2) introduce the four generating series $y_i(x) = \sum_{n \ge 0} a_{i,n} x^n$ where $a_{i,n} \in \{0, 1\}$ encodes for each $i$ a different condition on the parities of the number of occurrences of 1 and 2 in the ternary expansion of $n$. After deriving the linear Mahler system

$$Y(x) = \begin{pmatrix} 1 & x & 0 & x^2 \\ x & 1 & x^2 & 0 \\ 0 & x^2 & 1 & x \\ x^2 & 0 & x & 1 \end{pmatrix} Y(x^3) \quad \text{for} \quad Y(x) = \begin{pmatrix} y_1(x) \\ y_2(x) \\ y_3(x) \\ y_4(x) \end{pmatrix},$$

they prove that the $y_i(x)$ are linearly independent over $\overline{\mathbb{Q}}(x)$.

We return to these examples and discuss some additional ones in Section 10.2. Example 1.1(4) will also be our running example in §8; see Example 10.6.

## 1.2. First-order factors and hypergeometric solutions

The present article continues a line of work initiated in our earlier publication (CDDM 2018), to which we refer for more context. This work started when, back in December 2015, the second author asked about making effective a differential transcendence criterion introduced in (Dreyfus, Hardouin, and Roques 2018). This criterion boils down to determining the rational function solutions of a nonlinear Mahler equation analogous to the Riccati equation associated with a linear differential equation. Equivalently, this can be viewed as computing the first-order right-hand factors with rational function coefficients of the linear Mahler operator underlying the Mahler equation. We return to this motivation in §11 and postpone to that section a more detailed discussion of differential transcendence and differential algebraic independence criteria based on Mahler equations.

Another motivation for studying first-order factors is that, in the differential case, Beke's method (Markoff 1891; Bendixson 1892; Beke 1894)<fnref n="1" /> reduces the problem of factoring a differential equation into irreducible factors to that of finding the first-order right-hand factors of some auxiliary equations. One may expect that the method adapts to the Mahler case, so the present work paves the way to future factorization algorithms.

Our algorithms require solving linear Mahler equations in various domains. To this end, we build on results from (CDDM 2018). There, we worked with polynomials, rational functions, and series whose coefficients were in a computable subfield $\mathbb{K}$ of $\mathbb{C}$, so we continue with this assumption here. The following definitions and notation will be used throughout.

<footnote n="1">¹See (Bostan, Rivoal, and Salvy 2024) for more on the history of this method.</footnote>
**Definition 1.2.** For some fixed integer $b \ge 2$, called the radix, we denote by $M$ the Mahler operator with regard to $b$, which is defined as mapping any function $y(x)$ to $y(x^b)$. For $r \in \mathbb{N}$ and polynomials $\ell_i \in \mathbb{K}[x]$ with coefficients in a field $\mathbb{K}$ to be fixed by the context, we consider the linear Mahler equation

$$ \ell_r M^r y + \dots + \ell_1 My + \ell_0 y = 0, \qquad (L) $$

as well as the corresponding linear operator

$$ L = \ell_r M^r + \dots + \ell_1 M + \ell_0. $$

We assume $\ell_r \neq 0$ and call $r$ the order of the equation (L) and of the operator $L$. We write $d$ for the maximal degree $\max_k \deg \ell_k$ of the coefficients. We call

$$ \ell_r u M u \dots M^{r-1} u + \dots + \ell_2 u M u + \ell_1 u + \ell_0 = 0, \qquad (R) $$

the Riccati Mahler equation associated with (L).

**Hypothesis 1.3.** Throughout, we further assume that $\ell_0$ is also nonzero and that $L$ is primitive, that is, that the family of the $\ell_i$ has $\gcd 1$.

**Definition 1.4.** By a ramified rational function, we mean a rational function in some fractional power of $x$, so that rational functions are particular ramified rational functions.

Our goal is to find the right-hand factors $M - u$ of $L$ where $u$ is a ramified rational function. Equivalently, we want to compute the hypergeometric solutions of $L$, that is, the solutions $y$ of $L$ that satisfy a first-order linear Mahler equation with ramified rational coefficients. Informally, the link between both equations is the relation $u = My/y$. However, general solutions of (L) may live in a ring containing zero divisors, so that quotients are not always well defined (see the construction of $\mathfrak{D}$ in §3.1 and Definition 3.1).

**Example 1.5.** By searching for a linear dependency among the first coordinates of the iterates $Y, MY, M^2Y, \dots$ of the vector $Y$ appearing in Example 1.1(4), one obtains a linear Mahler operator of order 4 and degree 258 annihilating $y_1$. (See Example 10.2 in §10.2 for details of the method.) This operator turns out to annihilate all four $y_i$, as one can check by applying the method of Example 10.2 to each coordinate of the vector $Y$ in turn. The algorithms developed in the following sections reveal that its hypergeometric solutions are exactly the series

$$ \begin{gathered} a \prod_{k=0}^{\infty} (1 - x^{3^k} - x^{2 \cdot 3^k}), \qquad b \prod_{k=0}^{\infty} (1 + x^{3^k} - x^{2 \cdot 3^k}), \\ (c_0 + c_1 x) \prod_{k=0}^{\infty} (1 + x^{2 \cdot 3^k} + x^{4 \cdot 3^k}), \end{gathered} \qquad (1.1) $$

for arbitrary constants $a, b, c_0, c_1$. Comparing with the series expansions of $y_1, \dots, y_4$, we deduce that none of the $y_i$ is hypergeometric, let alone rational. By the classical dichotomy for Mahler functions (see §11), all the $y_i$ are therefore transcendental (i.e., not an algebraic function of $x$).
### 1.3. Contribution

This article provides algorithms to compute the rational solutions of Riccati Mahler equations, and more generally their ramified rational solutions. We develop two approaches.

The first one is adapted from the classical algorithm by Petkovšek (1992) for finding the hypergeometric solutions of a linear recurrence equation. Petkovšek's algorithm searches for first-order factors of difference operators in a special form called the Gosper–Petkovšek form. The existence of Gosper–Petkovšek forms relies on the fact that for any two nonzero polynomials $A$ and $B$, the set of integers $i$ such that $A(x)$ and $B(x+i)$ have a nontrivial common divisor is finite. As this is not true when the shift is replaced by the Mahler operator, we need to slightly depart from the classical definition of Gosper–Petkovšek forms. Doing so, we obtain a first complete algorithm for finding first-order factors of Mahler equations of arbitrary order. Note that Roques (2018) recently gave a slightly different adaptation of Petkovšek's algorithm to the case of Mahler equations of order two. Like Petkovšek's, these algorithms have to consider an exponential number of separate subproblems in the worst case. We discuss several ways of pruning the search space to mitigate this issue in practice.

Our second algorithm avoids the combinatorial search phase entirely, at the price of a worst-case exponential behavior of a different nature. It is based on a relaxation of the problem that can be solved using Hermite–Padé approximants<fnref n="2" />. The idea is to search for series solutions $y$ of the linear Mahler equation that make $My/y$ rational. Roughly speaking, we first compute a basis $(y_1, \dots, y_N)$ of series solutions, then search for linear combinations with polynomial coefficients of $y_1, \dots, y_N$, $My_1, \dots, My_N$ that vanish to a high approximation order $\sigma$, and finally isolate, among these relations, those corresponding to hypergeometric solutions by solving a polynomial system. Though we are not aware of any exact analogue of our algorithm in the literature, variants of the same basic idea have been used by several authors in the differential case (see below for references).

In order to state the algorithms and justify their correctness, it is useful to work in a ring containing “all” solutions of the linear equation (L), or at least all solutions needed in the discussion. Instead of appealing for this to the general Picard–Vessiot theory of linear difference equations, we introduce suitable conditions (Hypothesis 2.4) that suffice for our purposes and construct a ring satisfying them, whose elements (unlike those of Picard–Vessiot rings of Mahler equations) admit simple representations as formal series. In passing we define a suitable notion of Mahler-hypergeometric function and establish its basic properties.

We compare the performance of our two approaches based on an implementation<fnref n="3" /> in Maple and observe that the second algorithm turns out to be more

<footnote n="2">²We find some pleasant irony in our application of Hermite–Padé approximants to problems on Mahler operators, after Mahler himself has introduced similar approximants in his work (Mahler 1968).</footnote>

<footnote n="3">³Authored by Ph. Dumas. Available at https://mathexp.eu/dumas/dcfun/.</footnote>
efficient in practice on examples from the mathematical literature. Finally, we use this implementation in combination with criteria such as the one mentioned earlier to prove several differential algebraic independence results between series studied in the literature.

1.4. Related work

To the best of our knowledge, the problem we consider here was first discussed
in the doctoral dissertation of Dumas (1993, §3.6), which contains an incomplete
method for finding hypergeometric solutions of Mahler equations. Dumas'
method is somewhat reminiscent of Petkovšek's algorithm, but lacks a proper
notion of Gosper-Petkovšek form. Roques (2018, §6), as already mentioned,
describes a complete analogue of Petkovšek's algorithm for Mahler equations,
but restricts himself to equations of order two. We are not aware of any other
reference dealing specifically with the factorization of Mahler operators.

However, it is natural to try and adapt to Mahler equations algorithms that
apply to differential equations or to difference equations in terms of the shift
operator. In the differential case, factoring algorithms are a classical topic, going
back at least to Fabry's time (1885, §V); we refer to (Bostan, Rivoal, and Salvy
2024, §1.4) for a well-documented overview. The second of our algorithms, using
Hermite-Padé approximation, is related to methods known from the differential
case. Most similar to our approach is maybe an unpublished article by Bronstein
and Mulders (n.d. [1999?], §4) where they present a heuristic method for solving
Riccati differential equations based on Padé (not Hermite-Padé) approximants
of series solutions with indeterminate coefficients. The same idea appears in
works by Pflügel (1997, §2.5.2) and by van der Put and Singer (2003, §4.1),
though neither of these references discusses in detail how to deal with drops
in the degree of candidate solutions for special values of the parameters. The
core idea of enforcing the vanishing of high-order terms of series solutions of the
Riccati equation so as to reduce to the solution of polynomial equations in a
number of unknowns bounded by the order of the differential equation already
appears in Fabry's 1885 thesis.

In the shift case, the classical algorithm for finding hypergeometric solutions
is that of Petkovšek (1992). It is the direct inspiration for our first method.
Petkovšek's algorithm is itself based on Gosper's hypergeometric summation
algorithm (Gosper 1978), and was previously adapted to *q*-difference equations
in (Abramov and Petkovšek 1995; Abramov, Paule, and Petkovšek 1998). Van
Hoeij and Cluzeau (van Hoeij 1999; Cluzeau and Van Hoeij 2004) later proposed
alternative algorithms that are faster in practice; it seems likely that the ideas
would also be relevant in the Mahler case, but we do not explore this question here.

Another line of work aims at unified algorithms for linear functional operators
of various types based on the formalism of Ore polynomials (e.g., Bronstein
and Petkovšek 1993). Factoring algorithms, however, remain specific to each
individual type of equation. In addition, even methods that do apply to almost all
types of Ore operators sometimes fail in the Mahler case because the commutation
$Mx = x^bM$ does not preserve degrees with respect to $x$.
Turning now to the structure and computation of generalized series solutions of linear Mahler equations, Roques’s discussion of the local exponents of Mahler systems (2018, §4) forms the basis for our §3. Further developments (not used here) include recent work by Roques (2023), Faverjon and Roques (2024), and Faverjon and Poulet (2022).

1.5. *Outline*

The “structural” results about the solution spaces of Mahler equations come
first in the text. We first generalize our Mahlerian problem to the context
of difference rings, including the case of nonsurjective morphisms: in §2, we
derive a parametric description of the right-hand factors of a difference operator
(Theorem 2.9). To accommodate the Mahler operator in the previous generalized
framework, we then introduce in §3 an explicit difference ring $\mathcal{D}$ containing all
the solutions of the linear Mahler equations that are needed in the parametric
descriptions of the solutions $u$ to Riccati Mahler equations. In that section, we
also define classes of $F$-hypergeometric solutions for various difference fields $F$
and, for $F$ the field of Puiseux series, we partition the set of Puiseux series
solutions of the Riccati equation according to the coefficient of their valuation
term (Theorem 3.22). To describe the ramified rational solutions of the Riccati
equation, we then change $F$ to the field of ramified solutions in §4, where we
obtain a partition (Theorem 4.5) that is finer than the partition induced by
the partition in the Puiseux series on ramified solutions. In preparation for the
algorithms, the technical §5 presents bounds on the degree of numerators and
denominators of rational solutions to a Riccati Mahler equation (Proposition 5.2).

We continue with two algorithmic approaches.

We review Petkovšek's classical algorithm for the shift case and Roques's
analogue for order 2 in the Mahler case before developing our generalization in §6.
We first propose a brute-force algorithm (Algorithm 1), which we prove to be
correct (Theorem 6.5). Next, in §7 we propose several pruning criteria to improve
its exponential behavior, leading to an improved algorithm (Algorithm 3).

We then develop our relaxation method based on Hermite-Padé approxima-
tions: after studying in §8 the approximate syzygy module of a basis of truncated
solutions of (L) and its limit as precision goes to infinity, we obtain in §9 an
algorithm (Algorithm 4) that produces successive sets of parametrized candidates,
each containing all true solutions, until it stabilizes on true solutions only.

To conclude, we present examples of applications and timings in §10, where
we can observe that our relaxation method beats the other by combinatorial ex-
ploration in a number of natural examples. We finally apply our implementation
in order to prove properties of differential transcendence on examples in §11.

Part II on our generalized Petkovšek algorithm and Part III on our approxi-
mants-based algorithm are completely independent. They are both based on the
structural results of Part I. Part IV is mostly independent from the other parts
if one admits the existence of the various algorithms. Also, a few sections are
very independent from the rest of the text and might be skipped by readers with
a specific interest: §5 provides degree bounds on rational solutions that are to be
used as a black box in Part III; §7 presents technical algorithmic improvements to speed up our first plain generalization of Petkovšek's algorithm; §8.7 describes how approximate syzygies corresponding to nonsimilar (Definition 2.7) solutions interact; §9.2 sketches how the approach to the rational solving of the Riccati Mahler equation (R) adapts to the linear Mahler equation (L); §11 presents an application whose background in number theory is a bit away from our main theme, although it was our original motivation for the work.

## 1.6. Notation

The following is a collection of notation used throughout the text.

**Notation 1.6.** We will denote the field of Puiseux series over a field $\mathbb{L}$ by

$$ \mathbb{L}((x^{1/*})) = \bigcup_{q \ge 1} \mathbb{L}((x^{1/q})) $$

and the field of ramified rational functions over $\mathbb{L}$ by

$$ \mathbb{L}(x^{1/*}) = \bigcup_{q \ge 1} \mathbb{L}(x^{1/q}). $$

We use the notation $\mathbb{L}[x] \langle M \rangle$ to denote the algebra generated by $M$ over $\mathbb{L}[x]$ and satisfying the relation $Mf(x) = f(x^b)M$ for any $f \in \mathbb{L}[x]$. Similar definitions apply for other coefficient domains. All of the operators and functions $M$, val, ln, log take precedence over additions and products, which means $My/y = (My)/y$, val $ab = (val a)b$, etc. We write $Ly$ for the application of an operator $L$ to a function, or $L(x, M)y(x)$ if needed. We write $S_{\neq 0}$ for a given set $S$ containing $0$ to denote $S \setminus \{0\}$. The set $S$ will be the set $\mathbb{N}$ of natural numbers, a field, a vector space, a cone (see Definition 2.3), a Cartesian product, etc. A transpose is denoted by an exponent $T$. A tuple $z$ is identified with a row vector, its transpose $z^T$ with a column vector.

# Part I: Structural results

## 2. Structure of the solution set of the Riccati equation

The notions introduced in this section will be used when solving Riccati Mahler equations. As they do not depend on the specific choice of the Mahler operator, we write them in the broader generality of difference rings.

### 2.1. Basic notions of difference algebra

A difference ring is commonly defined as a pair $(D, \sigma)$ formed by a commutative ring with identity $D$ and an automorphism $\sigma$ of $D$; see, e.g., (Cohn 1965) or (van der Put and Singer 1997), and a difference field is a difference ring that is a field. An example is the field $\mathbb{K}((x^{1/*}))$ of Puiseux series, equipped with the Mahler operator $M$. However, we relax the definition to accept a map $\sigma$
that is only an injective endomorphism (cf. Wibmer 2013), since we intend to consider the restriction of the Mahler operator from $K(x)$ to itself, which is not an automorphism of the field of rational functions. When the context is clear, we write $D$ instead of $(D, \sigma)$.

Given a difference ring $(D, \sigma)$, a difference ring extension is a difference ring $(D', \sigma')$ such that $D'$ is an extension ring of $D$ and $\sigma'$ restricted to $D$ is equal to $\sigma$. In practice, we will always write again $\sigma$ for the extended $\sigma'$. The ring of constants of a difference ring $D$, denoted $D^\sigma$, is the subring of elements in $D$ fixed by the endomorphism $\sigma$. If $D$ is a difference field, the ring of constants is a field.

In what follows, every difference ring will be a difference ring extension of the field $K(x)$ of rational functions. For the needed level of generality, let us introduce

$$ \ell_r \sigma^r y + \dots + \ell_1 \sigma y + \ell_0 y = 0, \quad (L_\sigma) $$

$$ \ell_r u \sigma u \dots \sigma^{r-1} u + \dots + \ell_2 u \sigma u + \ell_1 u + \ell_0 = 0, \quad (R_\sigma) $$

where $\ell_0 \ell_r \neq 0$ and each $\ell_i$ is in $K(x)$. We will refer to these equations as, respectively, the linear difference equation $(L_\sigma)$ and the Riccati difference equation $(R_\sigma)$.

## 2.2. First-order factors and their solutions

We now explain the link between the linear and the Riccati equations: the solutions to the Riccati equation are the coefficients $u$ of the monic first-order right-hand factors $\sigma - u$ of the linear difference operator $L$.

**Lemma 2.1.** Given an element $u$ of a difference field extension $F$ of $K(x)$, the operator

$$ L = \ell_r(x)\sigma^r + \dots + \ell_0(x) \in K(x)\langle\sigma\rangle $$

associated with the linear difference equation $(L_\sigma)$ admits $\sigma - u$ as a first-order right-hand factor in $F\langle\sigma\rangle$ if and only if $u$ satisfies the Riccati difference equation $(R_\sigma)$.

*Proof.* The ring $F\langle\sigma\rangle$ is a skew Euclidean ring (Bronstein and Petkovšek 1993). By Euclidean division of $L$ on the right by $\sigma - u$ in the ring $F\langle\sigma\rangle$, we obtain $L = \tilde{L}(\sigma - u) + R$ where $R$ is exactly the left-hand side of $(R_\sigma)$. $\square$

We define hypergeometric elements in analogy with the classical difference case for shift operators.

**Definition 2.2.** Given a difference field extension $F$ of $K(x)$ and a difference ring extension $D$ of $F$, an element $y$ of $D$ is $F$-hypergeometric if there exists $u \in F$ satisfying $\sigma y = uy$.

Note that $y$ may be zero in the definition, but that we will focus on nonzero hypergeometric $y$ throughout.

The set of $F$-hypergeometric elements is generally not stable under addition, but it enjoys the structure of an $F$-cone, a notion that we introduce now.
<figure><img src="image_2.png" /><figcaption>Figure 1: The inclusion relations for a difference field extension $F$ of $K(x)$ and a ring extension $D$ of it satisfying Hypothesis 2.4. The two rings $F$ and $D$ share the same field of constants $L$. All inclusion but the dashed ones are difference ring inclusions; the dashed ones are only inclusions of $L$-vector spaces. The space $ker_D(L)$ is the space of the solutions of the operator $L$ underlying the equation $(L_\sigma)$.</figcaption></figure>

**Definition 2.3.** Given a field $F$, a nonempty set stable under multiplication by elements of $F$ will be called an $F$-cone.

The following hypothesis captures the notion of an extension of some difference field $F$ that contains “enough” $F$-hypergeometric solutions of the linear equation $(L_\sigma)$ to suitably describe the “full” solution set of the Riccati equation $(R_\sigma)$ in $F$. The hypothesis is illustrated in Figure 1.

**Hypothesis 2.4.** Let $F$ be a difference field extension of $K(x)$, with constant field $L$ containing $K$. A difference ring extension $D$ of $F$ is said to satisfy Hypothesis 2.4 if

1. the constant ring of $D$ is the field $L$,

2. for each nonzero $u \in F$, the equation $\sigma y = uy$ has nonzero solutions in $D$,

3. $D$ is simple, meaning that $\sigma$ is an automorphism and that the only ideals of $D$ stable under $\sigma$ are $(0)$ and $D$.

In contrast with the classical theory, we do not consider any universal Picard–Vessiot algebra (van der Put and Singer 1997, Prop. 1.33; Roques 2018, Theorem 35): by point 2, only equations of order 1 are known to have solutions in $D$. On the other hand, Hypothesis 2.4 ensures the following natural bound on the dimension of solution space.

**Lemma 2.5.** Let $D \supseteq F \supseteq K(x)$ be difference rings satisfying of Hypothesis 2.4 for a suitable constant field $L$ containing $K$.

Then, the extension $D$ is such that for any linear difference equation $(L_\sigma)$ of order $r$ with coefficients in $K(x)$, the vector space of solutions in $D$ of this equation has dimension at most $r$ over $L$.
*Proof.* By the assumption of simplicity (van der Put and Singer 1997, Def. 1.1; Singer 2016, Def. 4.1 and 4.4), the ring $D^\sigma$ of constants is a field (van der Put and Singer 1997, Lemma 1.7). By our assumption $\ell_0\ell_r \neq 0$ made in the introduction (Hypothesis 1.3), the companion matrix $A$ of $L$ is invertible over $K(x)$ and hence nonsingular over $D$. Suppose $y_1, \dots, y_m$ are $\mathbb{L}$-linearly independent solutions of $Ly = 0$, and, for $1 \le i \le m$, introduce $Y_i = (y_i, \sigma y_i, \dots, \sigma^{r-1}y_i)^T$, which satisfies $\sigma Y_i = AY_i$. By point 1 and (Singer 2016, Lemma 4.8), this family $\{Y_i\}_{i=1}^m$ of solutions of $\sigma Y = AY$, which is linearly independent over $\mathbb{L}$, is linearly independent over $D$ in the free $D$-module $D^r$. By the commutativity of $D$, the relation $m \le r$ holds, thus proving the lemma. $\square$

### 2.3. Parametrization of the solution set

Let $F$ be a fixed difference field extension of $K(x)$. Note that an extension $D$ of $F$ satisfying Hypothesis 2.4 is only a ring, so $\sigma y/y$ is not defined for all $y \in D$. However, restricting to nonzero $F$-hypergeometric $y$ makes $1/y$ well defined, as stated in the following lemma.

**Lemma 2.6.** Any nonzero $F$-hypergeometric element $y$ of some difference ring extension $D$ satisfying Hypothesis 2.4 over $F$ has an $F$-hypergeometric inverse in $D$.

*Proof.* For any nonzero $F$-hypergeometric $y \in D$, introduce $u \in F$ satisfying $\sigma y = uy$. The ideal generated by $y$ in $D$ contains $1$ by simplicity, so $y$ has an inverse $z$. Consequently, $\sigma z = 1/(\sigma y) = 1/(uz) = (1/u)z$, proving that $z$ is $F$-hypergeometric. $\square$

We now define similarity, a key concept in what follows.

**Definition 2.7.** Given a difference field extension $F$ of $K(x)$ and a difference ring extension $D$ of $F$, two elements $y_1$ and $y_2$ of $D$ are $F$-similar if there exists a nonzero $q \in F$ satisfying $y_1 = qy_2$.

Similarity is an equivalence relation, and $y_1$ and $y_2$ in the definition are either both nonzero or both zero. In what follows, we focus on nontrivial equivalence classes, that is, other than $\{0\}$, which therefore do not contain 0. We will denote such nontrivial similarity classes by $\mathfrak{H}_{\neq 0}$, reserving for the augmentation by 0 of these similarity classes the notation $\mathfrak{H} = \mathfrak{H}_{\neq 0} \cup \{0\}$. (Cf. Notation 1.6 in the introduction.)

**Lemma 2.8.** For any $F$-similarity class $\mathfrak{H}_{\neq 0}$ of $F$-hypergeometric elements, the set $\mathfrak{H}$ is an $\mathbb{L}$-vector space.

*Proof.* Because $\mathbb{L}$ is a subfield of $F$, $\mathfrak{H}$ is an $\mathbb{L}$-cone. We show that it is also stable under addition. Consider two elements $y_1$ and $y_2$ in $\mathfrak{H}_{\neq 0}$ such that $y_1 + y_2 \neq 0$, and $u$ and $q$ in $F_{\neq 0}$ such that $\sigma y_1 = uy_1$ and $y_2 = qy_1$. Then, as $q \neq -1 \neq \sigma q$, $\sigma(y_1 + y_2) = (1 + \sigma q)uy_1 = (1 + \sigma q)(1 + q)^{-1}u(y_1 + y_2)$, from which follows that $y_1 + y_2$ is $F$-hypergeometric and $F$-similar to $y_1$, thus in $\mathfrak{H}_{\neq 0}$. $\square$
We use the notation $\mathbb{P}(V)$ to denote the projectivization of a vector space $V$, leaving implicit the field over which $V$ is defined. We also write $\mathbb{P}^n(K)$ for the $K$-projective space $\mathbb{P}(K^{n+1})$ of dimension $n$ over a field $K$.

**Theorem 2.9.** Let $D$ be a difference ring extension of $F$ satisfying Hypothesis 2.4. For a fixed $(L_\sigma)$, the following holds:

1. The $L$-cone of $F$-hypergeometric solutions of $(L_\sigma)$ in $D$ naturally partitions into the class $\{0\}$ and nontrivial $F$-similarity classes $\mathfrak{S}_{\neq 0}$. Moreover, the $L$-cones $\mathfrak{S}$ are $L$-vector spaces in direct sum, and their dimensions add up to at most $r$.

2. The set of solutions in $F$ of $(R_\sigma)$ partitions into the images under the map $y \mapsto \sigma y/y$ of all nontrivial $F$-similarity classes of the set of $F$-hypergeometric solutions of $(L_\sigma)$.

3. Each nontrivial $F$-similarity class $\mathfrak{S}_{\neq 0} \subseteq D$ of the set of $F$-hypergeometric solutions of $(L_\sigma)$ induces a one-to-one parametrization of its image under $y \mapsto \sigma y/y$ by the $L$-projective space $\mathbb{P}(\mathfrak{S})$.

4. The sum of the dimensions of the parametrizing $L$-projective spaces $\mathbb{P}(\mathfrak{S})$ does not exceed the order $r$ of the equation $(L_\sigma)$ minus the number of nontrivial similarity classes $\mathfrak{S}_{\neq 0}$.

*Proof.* 1. For each $F$-similarity class $\mathfrak{S}_{\neq 0}$ of $F$-hypergeometric solutions of $L$, Lemma 2.8 implies that $\mathfrak{S}$ is an $L$-vector space. We now prove that all the $\mathfrak{S}$ are in direct sum. Given a relation $\lambda_0 y_0 + \cdots + \lambda_N y_N = 0$ with coefficients $\lambda_k \in L_{\neq 0}$ between nonzero $F$-hypergeometric solutions $y_k$ in distinct $F$-similarity classes, introduce $u_k \in F$ such that $\sigma y_k = u_k y_k$ for $0 \le k \le N$. Applying $\sigma - u_0$ yields

$$ \lambda_1(u_1 - u_0)y_1 + \cdots + \lambda_N(u_N - u_0)y_N = 0. \quad (2.1) $$

Because the $y_k$ were taken from distinct classes, each of the $y_0/y_k$, which exists by Lemma 2.6, cannot be in $F$, let alone be a constant. So, $\sigma(y_0/y_k) = (u_0/u_k) \times (y_0/y_k)$ is not $y_0/y_k$, and $u_0/u_k$ is not 1. Therefore, none of the coefficients $\lambda_k(u_k - u_0)$ for $1 \le k \le N$ is zero, and (2.1) is a shorter nontrivial relation. Iterating the process, we obtain that some $y_k$ is zero, a contradiction. As a consequence, the $L$-vector spaces $\mathfrak{S}$ are in direct sum. By the conclusion of Lemma 2.5, this direct sum is included in an $L$-space of dimension at most $r$, hence the dimension bound.

2. Given a solution $u \in F$ of $(R_\sigma)$, Lemma 2.1 implies that $\sigma - u$ is a right-hand factor of $L$. Because of Hypothesis 2.4, there exists $y \neq 0$ such that $\sigma y = uy$, and $y$ is invertible by Lemma 2.6. As a consequence, $y$ is also a solution of $(L_\sigma)$ and $u$ equals $\sigma y/y$. So the map $y \mapsto \sigma y/y$ is a (well-defined) surjection from the set of nonzero $F$-hypergeometric solutions of $(L_\sigma)$ to the set of solutions of $(R_\sigma)$ in $F$. If $u$ can be obtained by two solutions $y_1$ and $y_2$, that is, if $\sigma y_1/y_1 = \sigma y_2/y_2$, then $\sigma(y_1/y_2) = y_1/y_2$ is a constant in $D$, thus in $L$, so that, $y_1$ and $y_2$ are $F$-similar. As a consequence, the partitioning of the set of nonzero
F-hypergeometric solutions of (L_σ) translates into a partitioning of the set of solutions of (R_σ) in F.

3. By the proof of point 2, a nontrivial F-similarity class $\mathfrak{H}_{\neq 0}$ of the set of F-hypergeometric solutions of (L_σ) maps under $y \mapsto \sigma y/y$ to a subset U of the set of solutions of (R_σ) in F. Given $y_1$ and $y_2$ in the same class $\mathfrak{H}_{\neq 0}$, thus satisfying $y_1 = qy_2$ for some $q \in F$ (which needs to be nonzero), if both share the same image in U, that is, $\sigma y_1/y_1 = \sigma y_2/y_2$ is the same element of U, then the proof of 2 has shown that $q$ is a constant from $\mathbb{L}$. This implies that $y_1$ and $y_2$ are in the same $\mathbb{L}$-projective class, thus proving the result.

4. Point 1 implies that each similarity class $\mathfrak{H}$ is finite-dimensional. If it has $\mathbb{L}$-dimension $s$, its associated $\mathbb{L}$-projective space $\mathbb{P}(\mathfrak{H})$ has dimension $s-1$. The statement is next a direct consequence of points 1 and 3. □

Remark 2.10. Point 3 of Theorem 2.9 leads to a more explicit parametrization, owing to the finite $\mathbb{L}$-dimension of the space $\mathfrak{H}$. Fix $y_0$ in $\mathfrak{H}_{\neq 0}$ and let $u_0$ satisfy $\sigma y_0 = u_0 y_0$. Let $(q_1 y_0, \dots, q_s y_0)$ be some basis of $\mathfrak{H}$. With coordinates, the bijection of the theorem translates into the map

$$ (c_1 : \dots : c_s) \mapsto \frac{\sigma(c_1 q_1 + \dots + c_s q_s)}{c_1 q_1 + \dots + c_s q_s} u_0, $$

now a rational bijection from $\mathbb{P}^{s-1}(\mathbb{L})$ to its image.

Remark 2.11. The partition obtained by point 2 of the theorem can be interpreted as the partition induced on the solution set of (R_σ) by the equivalence relation on $F_{\neq 0}$ defined by: $u$ and $u'$ are F-equivalent if there exists $f \in F_{\neq 0}$ such that $u/u' = \sigma f/f$.

### 3. Generalized series solutions of the linear equation

In this section, we define a Mahler difference ring, denoted $\mathfrak{D}$ (Definition 3.1), of explicit formal expressions that is an extension satisfying Hypothesis 2.4 of both $\overline{\mathbb{K}}(x^{1/*})$ and $\overline{\mathbb{K}}((x^{1/*}))$. This will be a crucial ingredient to describe the Puiseux series solutions of a Riccati Mahler equation (Theorem 3.22). To this end, we describe $\mathbb{L}((x^{1/*}))$-hypergeometric elements (Proposition 3.18) so as to apply Theorem 2.9 for a given field $\mathbb{L}$. The properties of $\mathfrak{D}$ are again used in §4 to describe the ramified rational solutions of the Riccati equation (Theorem 4.5).

#### 3.1. The difference ring $\mathfrak{D}$

Slightly adapting (Roques 2018, §5.2), we introduce the $\overline{\mathbb{K}}((x^{1/*}))$-algebra with basis $(e_\lambda)_{\lambda \in \overline{\mathbb{K}}_{\neq 0}}$ for elements $e_\lambda$ satisfying the algebraic relations $e_\lambda e_{\lambda'} = e_{\lambda \lambda'}$ and $e_1 = 1$. This ring is not a domain, as shown by the product $(e_\lambda + e_{-\lambda}) \times (e_\lambda - e_{-\lambda}) = 0$. Next, we equip it with a structure of difference ring by enforcing $Me_\lambda = \lambda e_\lambda$. To support the intuition, we henceforth denote $e_\lambda$ by $(\ln x)^{\log_b \lambda}$ for $\lambda \in \overline{\mathbb{K}}_{\neq 0}$. This mere notation bears no analytic meaning, but it is reminiscent of the relation $\ln(x^b) = b \ln(x)$ in analysis, as well as of the linear independence over the field of meromorphic functions of the family of functions $((\ln x)^{\log_b \lambda})_{\lambda \in \mathbb{C}_{\neq 0}}$.
**Definition 3.1.** We define a difference ring by

$$ \mathcal{D} = \bigoplus_{\lambda \in \bar{\mathbb{K}} \setminus \{0\}} (\ln x)^{\log_b \lambda} \bar{\mathbb{K}}((x^{1/\lambda})) \quad (3.1) $$

where $(\ln x)^{\log_b \lambda}$ is an eigenvector of $M$ with respect to the eigenvalue $\lambda$.

Given a subfield $\mathbb{L}$ of $\bar{\mathbb{K}}$, we also define a difference subring $\mathcal{D}_{\mathbb{L}}$ of $\mathcal{D}$ by

$$ \mathcal{D}_{\mathbb{L}} = \bigoplus_{\lambda \in \mathbb{L} \setminus \{0\}} (\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/\lambda})), \quad (3.2) $$

so that $\mathcal{D}$ is the special case $\mathcal{D}_{\bar{\mathbb{K}}}$.

For an operator of $\bar{\mathbb{K}}(x)\langle M \rangle$,

$$ L = L(x, M) = \sum_{k=0}^{r} \ell_k(x) M^k, \quad (3.3) $$

and an element of $\mathcal{D}$,

$$ y = \sum_{\lambda \in S} (\ln x)^{\log_b \lambda} p_{\lambda}(x), \quad (3.4) $$

where the set $S$ is a finite subset of $\bar{\mathbb{K}}$ and the $p_{\lambda}(x)$ are Puiseux series, the action of $L$ on $y$ has a special structure given by

$$ L(x, M) y = \sum_{\lambda \in S} (\ln x)^{\log_b \lambda} L(x, \lambda M) p_{\lambda}(x). \quad (3.5) $$

This formula decomposes $(\mathbb{L})$ into equations on the $p_{\lambda}(x)$ as follows.

**Lemma 3.2.** For an operator $L \in \bar{\mathbb{K}}(x)\langle M \rangle$ as in (3.3) and an element $y$ of $\mathcal{D}$ as in (3.4), the equality $L(x, M)y = 0$ is equivalent to the finite set of equalities

$$ L(x, \lambda M) p_{\lambda}(x) = 0, \quad \lambda \in S. \quad (3.6) $$

*Proof.* This follows directly from (3.5), by the uniqueness of coefficients in $\mathcal{D}$. $\square$

The following lemma can be deduced from Roques (2018, Theorem 35) in the case $\mathbb{L} = \bar{\mathbb{Q}}$.

**Lemma 3.3.** The field of constants of the ring $\mathcal{D}_{\mathbb{L}}$ is $\mathbb{L}$.

*Proof.* By Lemma 3.2, the equation $My = y$ is equivalent to $\lambda Mp_{\lambda} = p_{\lambda}$ for $\lambda \in \bar{\mathbb{K}}$. For each $\lambda$, if $p_{\lambda}$ is nonzero, its valuation must be 0 and its leading coefficient $c$ must satisfy $\lambda c = c$, implying $\lambda = 1$. This means that $y$ is in $\mathbb{L}$. We conclude that the constant field of $\mathcal{D}_{\mathbb{L}}$ is $\mathbb{L}$ itself. $\square$
## 3.2. Bounds

The goal of this section is to bound the valuation and ramification order of the Puiseux series solutions of (L) in $\mathfrak{D}$. We need to introduce notions related to a lower Newton polygon attached to L, which the reader should not confuse with parallel notions attached to the upper Newton polygon that will be defined in §6.3 (Definition 6.2).

**Definition 3.4.** To the polynomial

$$L = \sum_{k=0}^{r} \ell_k(x) M^k = \sum_{k,j} \ell_{k,j} x^j M^k$$

we associate the set of points $(b^k, j) \in \mathbb{R}^2$ such that $\ell_{k,j} \neq 0$. The lower convex hull of this set is called the lower Newton polygon (of L).

The following lemma is borrowed from (CDDM 2018, Lemma 2.2).

**Lemma 3.5.** The valuation of any Puiseux series solution of (L) is the opposite of the slope of an edge of the lower Newton polygon of L. Moreover the coefficients $\ell_{k,j}$ corresponding to all points $(b^k, j)$ lying on the edge must add up to zero.

For a solution $y$ of (L), that is, satisfying $L(x, M)y = 0$, to be of the form $(\ln x)^{\log_b \lambda} p(x)$ with $p(x)$ a Puiseux series, the series $p(x)$ must by the formula (3.5) be a solution of $L(x, \lambda M)p = 0$. The change from $L(x, M)$ to $L(x, \lambda M)$ does not modify the Newton polygon, but the coefficients at abscissa $b^k$ are multiplied by $\lambda^k$. So, by Lemma 3.5, the only $\lambda$ for which a nonzero solution of the form $(\ln x)^{\log_b \lambda} p(x)$ may exist are the (necessarily nonzero) roots of the characteristic polynomials defined by (3.7) in the following definition.

**Definition 3.6.** The edges of the lower Newton polygon of L are indexed so that their slopes form an increasing sequence. Given the abscissae $b^{k_1}$ and $b^{k_2}$, with $k_1 \le k_2$, of the jth edge, the polynomial $\chi_j$ defined by

$$\chi_j(X) := \sum_{k=k_1}^{k_2} \ell_{k,j} X^{k-k_1} \qquad (3.7)$$

is called the characteristic polynomial of the jth edge. We define $\Lambda$ as the union of the roots of these polynomials:

$$\Lambda = \{\lambda \in \bar{\mathbb{K}} : \exists j, \chi_j(\lambda) = 0\}. \qquad (3.8)$$

Given $\lambda$ in $\Lambda$, the jth edge will be called $\lambda$-admissible (for L) if $\chi_j(\lambda) = 0$.

The present 1-admissible edges were simply called admissible in (CDDM 2018). Lemma 3.5 extends slightly to accommodate a logarithmic part, leading to the following result. The omitted proof parallels that of (CDDM 2018, Lemma 2.2).
**Lemma 3.7.** The valuation of the Puiseux series $p(x)$ in any nonzero solution of (L) of the form $(\ln x)^{\log_b \lambda} p(x)$ is the opposite of the slope of a $\lambda$-admissible edge of the lower Newton polygon of L.

As already noted for $\lambda = 1$ in (CDDM 2018), $\lambda$-admissibility is only a necessary condition for the existence of a series solution with logarithmic part $(\ln x)^{\log_b \lambda}$.

**Proposition 3.8.** For an operator $L \in \bar{\mathbb{K}}(x)\langle M \rangle$ as in (3.3) and an element $y = (\ln x)^{\log_b \lambda} p(x)$ of $\mathfrak{D}$ given by $\lambda \in \bar{\mathbb{K}}_{\neq 0}$ and $p \in \bar{\mathbb{K}}((x^{1/*}))$, the equality $L(x, M)y = 0$ is equivalent to the equality

$$L(x, \lambda M) p(x) = 0. \qquad (3.9)$$

Moreover, (3.9) can have a nonzero solution only if $\lambda$ is the root of some characteristic polynomial of L.

*Proof.* The stated equivalence is a mere specialization of Lemma 3.2. The necessary condition on $\lambda$ has been proved above when introducing characteristic polynomials. $\square$

In order to discuss ramifications, we import (CDDM 2018, Prop. 2.19), which deals with 1-admissible edges. In doing so, we slightly reformulate the result so as to describe solutions over an algebraic extension $\mathbb{L}$ of $\mathbb{K}$.

**Lemma 3.9.** Any Puiseux series solution of (L) is an element of $\bar{\mathbb{K}}((x^{1/q_1}))$, where $q_1$ is the lcm of the denominators of the slopes of those 1-admissible edges such that said denominators are coprime with b. Furthermore, a $\bar{\mathbb{K}}$-vector basis of the solutions can be written in $\mathbb{K}((x^{1/q_1}))$. Given an intermediate field $\mathbb{L}$, such a $\bar{\mathbb{K}}$-vector basis is also an $\mathbb{L}$-vector basis for the $\mathbb{L}$-space of solutions in $\mathbb{L}((x^{1/*}))$.

*Proof.* The first assertion of the statement is exactly the proposition we referred to, using $\bar{\mathbb{K}}$ as a base field. The second and third assertions are implicit in (CDDM 2018), which finds solutions by solving linear systems over $\mathbb{K}$. $\square$

**Definition 3.10.** More generally, we denote by $q_\lambda$ the lcm of the denominators coprime with b of the slopes of the $\lambda$-admissible edges of the lower Newton polygon of L, that is

$$q_{\lambda} = \operatorname{lcm} \{q \in \mathbb{N} : \exists j, p, \operatorname{gcd}(b, q) = \operatorname{gcd}(p, q) = 1, \chi_j(\lambda) = 0, \text{ and the jth edge has slope } p/q\}. \quad (3.10)$$

We call $q_\lambda$ the ramification bound with respect to $\lambda$.

We can then prove the following result on solutions in $\mathfrak{D}_{\mathbb{L}}$ (Definition 3.1).

**Proposition 3.11.** Given an intermediate field $\mathbb{L}$, the solution space of (L) in $\mathfrak{D}_{\mathbb{L}}$ is included in

$$\bigoplus_{\lambda \in \Lambda \cap \mathbb{L}} (\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/q_\lambda})).$$
Furthermore, an $\mathbb{L}$-basis of this space can be written in

$$ \bigcap_{\lambda \in \Lambda \cap \mathbb{L}} (\ln x)^{\log_b \lambda} \mathbb{K}[\lambda]((x^{1/q_\lambda})). $$

*Proof.* For a nonzero solution $y$ of $L$ in $\mathcal{D}_{\mathbb{L}}$, Lemma 3.2 splits the relation $Ly = 0$ into the equivalent finite set of equations (3.6). This proves the first part, where the direct sum of the solution space is induced by the direct sum in the expression (3.2) of $\mathcal{D}_{\mathbb{L}}$. Next, according to Lemma 3.5, the equation (3.6) for a given $\lambda \in \mathbb{L}$ has a nonzero solution only if there exists a $\lambda$-admissible edge, that is, only if $\lambda$ is an element of $\Lambda$. Moreover, the change from $L(x, M)$ to $L(x, \lambda M)$ modifies the linear system that is used in the proof of Lemma 3.9 to determine the solutions by moving its coefficients to the algebraic extension $\mathbb{K}[\lambda]$. As a result, a basis of the Puiseux-series solutions $p_\lambda$ can be found with coefficients in $\mathbb{K}[\lambda]$ and ramification order (bounded by) $q_\lambda$. The result is proved. $\square$

The following corollary is already proved as Lemma 2.5: we obtain here an alternative proof by considering properties of the lower Newton polygon.

**Corollary 3.12.** Given an intermediate field $\mathbb{L}$, the solution set in $\mathcal{D}_{\mathbb{L}}$ of the linear equation (L) is an $\mathbb{L}$-vector space whose dimension over the field $\mathbb{L}$ is at most the order $r$ of the equation.

*Proof.* Proposition 3.11 provides us with an $\mathbb{L}$-basis of the solutions of (L) in $\mathcal{D}_{\mathbb{L}}$, which can be chosen such that for each $\lambda$ the basis elements of the form $(\ln x)^{\log_b \lambda} p(x)$ involve series $p(x)$ with distinct valuations. By Lemma 3.7, these valuations are opposites of slopes of distinct edges of the Newton polygon whose characteristic polynomials $\chi_j$ all have $\lambda$ as a root. Therefore, the total number of elements of the basis is bounded by $\sum_j \deg \chi_j \le r$. $\square$

The order $r$ is not always reached in Corollary 3.12, as shown by the following two examples.

**Example 3.13.** For any radix $b \ge 2$, the operator $L = M^2 - 2M + 1$ admits the series solution 1, and no nonconstant solution. This reflects the fact that the nonhomogeneous equation $My - y = 1$ has no solution in $\mathcal{D}$, while any solution of this equation needs to be killed by $L$.

**Example 3.14.** Similarly, when $b \ge 2$, $x^{b-1}M^2 - (1 + x^{b-1})M + 1$ admits the series solution 1, and otherwise would require a Hahn series $x^{-1/b} + x^{-1/b^2} + x^{-1/b^3} + \dots$ to get a second dimension of solutions, the latter solving $My - y = 1/x$. See (Faverjon and Roques 2024) for a definition and developments.

### 3.3. Puiseux-hypergeometric solutions

In the present section, we focus on the field $F = \mathbb{L}((x^{1/*}))$ to obtain the Puiseux series with coefficients in $\mathbb{L}$ that are solutions of the Riccati equation. We describe the $\mathbb{L}((x^{1/*}))$-hypergeometric elements in $\mathcal{D}_{\mathbb{L}}$ and deduce that $\mathcal{D}_{\mathbb{L}}$ satisfies Hypothesis 2.4 over $\mathbb{L}((x^{1/*}))$.
**Lemma 3.15.** Fix $q \in \mathbb{N}_{\neq 0}$ and $u \in \mathbb{L}((x^{1/q}))$ of valuation $\alpha$ and with leading coefficient $\lambda$. If $\lambda \neq 1$, then the equation $Mv = uv$ admits no nonzero solution in $\mathbb{L}((x^{1/*}))$. Otherwise, it admits a unique series solution $v$ in $x^{\alpha/(b-1)}\mathbb{L}[[x^{1/q}]]$ of valuation $\alpha/(b-1)$ and with leading coefficient 1, and the solution space of this equation in $\mathbb{L}((x^{1/*}))$ is exactly the one-dimensional space $\mathbb{L}v$.

*Proof.* The Puiseux series $u$ can be written in the form $u = \lambda x^{\alpha} \hat{u}$ for some $\hat{u} = \sum_{k \in \mathbb{N}} \hat{u}_k x^{k/q} \in \mathbb{L}[[x^{1/q}]]$ such that $\hat{u}_0 = 1$. Assume that the equation $Mv = uv$ admits a nonzero solution $v \in \mathbb{L}((x^{1/*}))$ with valuation $\beta$ and leading coefficient $v_0 \neq 0$. Taking leading terms on both sides yields $v_0 x^{b\beta} = \lambda v_0 x^{\alpha+\beta}$, which implies $\lambda = 1$. This proves the absence of nonzero solutions if $\lambda \neq 1$ and the announced valuation of $v$ if $\lambda = 1$.

We continue with the latter case. Write $v$ in the form $x^{\alpha/(b-1)}\hat{v}$ for a series $\hat{v} = \sum_{k \in \mathbb{N}} \hat{v}_k x^{k/q}$ and $\hat{v}_0 = v_0 \neq 0$, the equation is equivalent to the relation

$$ \sum_{k \in \mathbb{N}} \hat{v}_k x^{bk/q} = \sum_{n \in \mathbb{N}} \sum_{i \in \mathbb{N}} \hat{u}_{n-i} \hat{v}_i x^{n/q}, $$

which in turn is equivalent to the recurrence formula

$$ \hat{v}_n = \hat{v}_{n/b} - \sum_{i=0}^{n-1} \hat{u}_{n-i} \hat{v}_i $$

for $n \in \mathbb{N}_{\neq 0}$, where $\hat{v}_{n/b}$ is understood to be zero if $n/b$ is not in $\mathbb{N}$. The sequence $(\hat{v}_n)_{n \in \mathbb{N}}$, and therefore the solution $v$, is fully determined by the choice of $\hat{v}_0 = 1$.

Next, for any Puiseux series solution $w$ of the equation $Mw = uw$, the quotient $w/v$ is in the constant field of $\mathbb{L}((x^{1/*}))$, that is, it lies in $\mathbb{L}$. $\square$

The following result is to be found as (Roques 2018, Theorem 35, top of page 342), where $\mathfrak{D}$ is named $B$ and is defined in (Roques 2018, Lemma 36). We reproduce it here for the readers' convenience, while removing the unnecessary hypothesis of an algebraically closed field $\mathbb{L}$.

**Lemma 3.16.** $\mathfrak{D}_{\mathbb{L}}$ is a simple difference ring.

*Proof.* Let $J$ be a nonzero difference ideal of $\mathfrak{D}_{\mathbb{L}}$. Let $y = \sum_{\lambda \in \mathbb{L}_{\neq 0}} (\ln x)^{\log_b \lambda} p_\lambda$ for $p_\lambda \in \mathbb{L}((x^{1/*}))$ be a nonzero element of $J$ such that the cardinality of the support of $(p_\lambda)_{\lambda \in \mathbb{L}_{\neq 0}}$ is minimal. Let $\lambda_0$ be such that $p_{\lambda_0} \neq 0$ and replace $y$ with $y/p_{\lambda_0}$, so that we assume that $p_{\lambda_0} = 1$. Because $p_{\lambda_0}$ appears with zero coefficient in $y - My \in J$ and by the minimality assumption, $y - My = 0$. Consequently, for all $\lambda \in \mathbb{L}_{\neq 0}$, $p_\lambda = \lambda Mp_\lambda$. For $\lambda \neq 1$, we conclude $p_\lambda = 0$ by Lemma 3.15. Consequently, $y$ reduces to $p_1$, which is invertible, thus proving $J = \mathfrak{D}_{\mathbb{L}}$. $\square$

**Corollary 3.17.** $\mathfrak{D}_{\mathbb{L}}$ is an extension satisfying Hypothesis 2.4 over any intermediate difference field $F$ between $\mathbb{L}(x)$ and $\mathbb{L}((x^{1/*}))$, so that Theorem 2.9 applies to $\mathfrak{D}_{\mathbb{L}}$, viewed as an extension of $F$.

*Proof.* This follows directly from Lemma 3.3, Lemma 3.15, and Lemma 3.16. $\square$
We will use $F = \mathbb{L}((x^{1/*}))$ in Theorem 3.22 below, and $\mathbb{L}(x^{1/*})$ in Theorem 4.5.

**Proposition 3.18.** The $\mathbb{L}((x^{1/*}))$-hypergeometric elements in $\mathcal{D}_{\mathbb{L}}$ are its elements $y = (\ln x)^{\log_b \lambda} p(x)$ where $\lambda$ is in $\mathbb{L}_{\neq 0}$ and $p(x)$ is a Puiseux series in $\mathbb{L}((x^{1/*}))$. Two such nonzero elements are $\mathbb{L}((x^{1/*}))$-similar if and only if they share the same $\lambda$.

*Proof.* Let $y$ be a nonzero $\mathbb{L}((x^{1/*}))$-hypergeometric element in $\mathcal{D}_{\mathbb{L}}$, together with the corresponding Puiseux series $u$ satisfying the equation $My = uy$. The series $u$ lies in a subfield $\mathbb{L}((x^{1/q}))$ with $q$ a positive integer. Let its leading term be denoted $\lambda_0 x^{m/q}$, for suitable $\lambda_0 \in \mathbb{L}$ and $m \in \mathbb{Z}$. By Lemma 3.2, writing $y$ as in (3.4) for a finite subset $S$ of $\mathbb{L}_{\neq 0}$ and some $p_\lambda(x)$ in $\mathbb{L}((x^{1/*}))$ leads to the equations $\lambda M p_\lambda = u p_\lambda$ for $\lambda \in S$. Any such equation can have a nonzero solution only if $\lambda$ is equal to $\lambda_0$, forcing $y$ to be of the form $y = (\ln x)^{\log_b \lambda_0} p_{\lambda_0}(x)$ where $p_{\lambda_0}(x) \in \mathbb{L}((x^{1/*}))$ is nonzero. Conversely, any $y$ of that form is $\mathbb{L}((x^{1/*}))$-hypergeometric in $\mathcal{D}_{\mathbb{L}}$.

Now suppose that $y_1 = (\ln x)^{\log_b \lambda_1} p_1(x)$ and $y_2 = (\ln x)^{\log_b \lambda_2} p_2(x)$ are $\mathbb{L}((x^{1/*}))$-similar and nonzero. The similarity implies the linear dependence of $(\ln x)^{\log_b \lambda_1}$ and $(\ln x)^{\log_b \lambda_2}$ over $\mathbb{L}((x^{1/*}))$, and therefore the equality $\lambda_1 = \lambda_2$. Conversely two nonzero hypergeometric solutions of $(\mathbb{L})$ in $\mathcal{D}_{\mathbb{L}}$ with the same $\lambda \in \mathbb{L}$ are similar by definition. $\square$

*Remark 3.19.* The same approach, supplemented with Lemma 3.15, shows more precisely that the $\mathbb{L}((x^{1/q}))$-hypergeometric elements in $\mathcal{D}_{\mathbb{L}}$ are of the form

$$y = (\ln x)^{\log_b \lambda} x^{m/(q(b-1))} s(x^{1/q}) \qquad (3.11)$$

with $\lambda$ in $\mathbb{L}_{\neq 0}$, $q$ a positive integer, $m$ an integer, and $s(x)$ a series of valuation 0 in $\mathbb{L}[[x]]$. Observe (e.g., by making $m = 1$) that $My/y$ may have a lower ramification order than $y$. Consequently, obtaining all Laurent power series solutions to (R) in the form $u = My/y$ generally requires to consider ramified hypergeometric solutions $y$ of (L). A simple example is provided by $L = M - u$ for $b = 5$ and $u = x(1 + x^3)$, with solution $y = x^{1/4} \prod_{k \geq 0} (1 + x^{3 \cdot 5^k})^{-1}$.

*Remark 3.20.* We could have written the series $s$ in Proposition 3.18 as an infinite product. Indeed, this representation will prove useful later: in Proposition 4.1 where we obtain an infinite product of rational functions; in Proposition 4.3 where we provide a characterization of the similarity classes.

**Definition 3.21.** Let $\Lambda'$ denote the finite set of $\lambda \in \Lambda$ (see (3.8)) such that there exists a nonzero solution of $(\mathbb{L})$ in $(\ln x)^{\log_b \lambda} \mathbb{K}((x^{1/*}))$. For an intermediate field $\mathbb{L}$ and $\lambda \in \Lambda' \cap \mathbb{L}$, define $\mathfrak{H}_{\mathbb{L},\lambda}$ to be the set of hypergeometric solutions in $(\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/*}))$. Let $\mathfrak{R}_{\mathbb{L},\lambda} \subseteq \mathbb{L}((x^{1/*}))$ denote the image of $(\mathfrak{H}_{\mathbb{L},\lambda})_{\neq 0}$ under the map $y \mapsto My/y$.

The elements of $\mathfrak{H}_{\mathbb{L},\lambda}$ are necessarily in $(\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/q_\lambda}))$ by Proposition 3.11. With this notation, Theorem 2.9 specializes as follows.
**Theorem 3.22.** The solution set $\mathfrak{R}_{\mathbb{L}((x^{1/\lambda}))}$ of the Riccati equation (R) in the ramified rational function field $\mathbb{L}((x^{1/\lambda}))$ is the disjoint union

$$ \mathfrak{R}_{\mathbb{L}((x^{1/\lambda}))} = \coprod_{\lambda \in \mathbb{L} \cap \Lambda'} \mathfrak{R}_{\mathbb{L}, \lambda} \qquad (3.12) $$

indexed by the finite set $\mathbb{L} \cap \Lambda'$. Each $\mathfrak{R}_{\mathbb{L},\lambda}$ is a set of series from $\mathbb{L}((x^{1/q_{\lambda}}))$ with leading coefficient $\lambda$ and is one-to-one rationally parametrized by the $\mathbb{L}$-projective space $\mathbb{P}(\mathfrak{H}_{\mathbb{L},\lambda})$. The corresponding parametrization is obtained by restricting the map $y \mapsto My/y$ from $(\mathfrak{H}_{\mathbb{L},\lambda})_{\neq 0}$ to its image in $\mathbb{L}((x^{1/\lambda}))$ and projectivizing its source. Moreover the dimensions of the $\mathbb{L}$-projective spaces add up to a number that is at most the order $r$ of the linear equation (L) minus the cardinality of $\mathbb{L} \cap \Lambda'$.

*Proof.* By Corollary 3.17, the ring $\mathfrak{D}_{\mathbb{L}}$ satisfies Hypothesis 2.4 as an extension of $\mathbb{L}((x^{1/\lambda}))$, so Theorem 2.9 applies to $D = \mathfrak{D}_{\mathbb{L}}$ and $F = \mathbb{L}((x^{1/\lambda}))$. Points 1 to 3 of that theorem justify the partitioning according to (3.12) as well as the form of the parametrization. Moreover, from the inclusion $\mathfrak{H}_{\mathbb{L},\lambda} \subseteq (\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/q_{\lambda}}))$ follows the inclusion $\mathfrak{R}_{\mathbb{L},\lambda} \subseteq \mathbb{L}((x^{1/q_{\lambda}}))$. Point 4 proves the bound on dimensions. $\square$

# 4. Ramified rational solutions to the Riccati equation

We now study the $F$-hypergeometric solutions of (L) for $F = \mathbb{L}(x^{1/\lambda})$, or, equivalently, the ramified rational solutions of (R) with coefficients in some intermediate field $\mathbb{L}$.

## 4.1. Hypergeometric elements

We begin by characterizing the $\mathbb{L}(x^{1/\lambda})$-hypergeometric elements of $\mathfrak{D}_{\mathbb{L}}$ (Definitions 2.2 and 3.1).

**Proposition 4.1.** The nonzero $\mathbb{L}(x^{1/\lambda})$-hypergeometric elements in $\mathfrak{D}_{\mathbb{L}}$ are the elements

$$ y = c (\ln x)^{\log_b \lambda} x^{m/q} f(x^{1/q}) \qquad (4.1) $$

with $c$ and $\lambda$ in $\mathbb{L}_{\neq 0}$, $m$ an integer, $q$ a positive integer, and $f(x)$ a formal power series that is expressible as an infinite product

$$ f(x) = \prod_{k \geq 0} 1/g(x^{b^k}) \in \mathbb{L}[[x]] \qquad (4.2) $$

for some $g(x) \in \mathbb{L}(x)$ satisfying $g(0) = 1$.

*Proof.* Any nonzero $\mathbb{L}(x^{1/\lambda})$-hypergeometric element $y$ in $\mathfrak{D}_{\mathbb{L}}$ is $\mathbb{L}((x^{1/\lambda}))$-hypergeometric. By Proposition 3.18, it can be written $y = (\ln x)^{\log_b \lambda} x^{m/q} s(x^{1/q})$ with $\lambda$ in $\mathbb{L}_{\neq 0}$, $q$ a positive integer, $m$ an integer, and $s(x)$ in $\mathbb{L}[[x]]$ of valuation 0. One then has $My/y = \lambda x^{(b-1)m/q} s(x^{b/q})/s(x^{1/q}) \in \mathbb{L}((x^{1/\lambda})) \cap \mathbb{L}(x^{1/\lambda})$, so $g(x) := s(x^b)/s(x)$ is a rational function. The infinite product (4.2) for this value of $g$ converges to $s$, hence $y$ is of the form (4.1) with $f = s$. $\square$
Remark 4.2. If the hypergeometric element $y$ is such that $My/y$ is in $\mathbb{L}((x^{1/q'}))$, then $q = q'(b-1)$ can be used in formula (4.1), as a consequence of Remark 3.19.

## 4.2. Similarity classes

The set of $\mathbb{L}(x^{1/*})$-hypergeometric solutions of (L) decomposes into $\mathbb{L}(x^{1/*})$-similarity classes that refine the $\mathbb{L}((x^{1/*}))$-similarity classes $\mathfrak{H}_{\mathbb{L},\lambda}$ described in §3 (Definition 3.21). The following proposition gives a normal form for these new similarity classes. Recall from Proposition 3.11 that the ramification bound $q_{\lambda}$ defined by (3.10) is a bound on the ramification orders of solutions of (L) in $(\ln x)^{\log_{b} \lambda} \mathbb{L}((x^{1/*}))$, and recall as well the definition of $\Lambda'$ in Definition 3.21.

**Proposition 4.3.** Any nontrivial $\mathbb{L}(x^{1/*})$-similarity class $\mathfrak{H}_{\neq 0}$ of $\mathbb{L}(x^{1/*})$-hypergeometric solutions of (L) uniquely determines a tuple $(\lambda, g, (p_j)_{1 \le j \le s})$ made of:

* $\lambda$ in $\mathbb{L} \cap \Lambda'$,

* an integer $s > 0$ and Laurent polynomials $p_1, \dots, p_s$ of $\mathbb{L}[x, x^{-1}]$,

* a rational function $g$ in $\mathbb{L}(x)$ satisfying $g(0) = 1$,

and such that:

1. the $s$ generalized series

$$ (\ln x)^{\log_b \lambda} \times p_j(x^{1/q_\lambda}) \times \prod_{k \ge 0} \frac{1}{g(x^{b^k/q_\lambda})}, \quad 1 \le j \le s, \quad (4.3) $$

form a basis of the $\mathbb{L}$-vector space $\mathfrak{H}$,

2. the family $(p_1, \dots, p_s)$ is in reduced echelon form w.r.t. ascending degree (meaning in particular that the coefficient of minimal degree of $p_j$ is 1),

3. the elements $p_1, \dots, p_s$ are coprime in $\mathbb{L}[x, x^{-1}]$ (in other words, the $p_j/x^{\text{val}p_j}$ are coprime in $\mathbb{L}[x]$).

**Proof.** The $\mathbb{L}(x^{1/*})$-similarity class $\mathfrak{H}_{\neq 0}$ of $\mathbb{L}(x^{1/*})$-hypergeometric solutions of (L) is contained in an $\mathbb{L}((x^{1/*}))$-similarity class of $\mathbb{L}((x^{1/*}))$-hypergeometric elements of $\mathfrak{D}_{\mathbb{L}}$, which, by Proposition 3.18, is of the form $(\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/*}))_{\neq 0}$ for some $\lambda \in \mathbb{L}_{\neq 0}$. Since, additionally, $\mathfrak{H}$ is contained in the solution space of (L), Proposition 3.11 implies $\mathfrak{H} \subseteq (\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/q_\lambda}))$. Additionally, $\mathfrak{H}$ has finite dimension over $\mathbb{L}$ by Theorem 2.9.

For any element $y \in \mathfrak{H}_{\neq 0}$, use Proposition 4.1 to write

$$ y = c (\ln x)^{\log_b \lambda} x^{m/q} f(x^{1/q}), \quad f(x) = \prod_{k \ge 0} \frac{1}{g(x^{b^k})}, \quad (4.4) $$

with $c \in \mathbb{L}$, $g(x) \in \mathbb{L}(x)$, and $g(0) = 1$. We contend that we can without loss of generality make $q$ equal to the ramification bound $q_{\lambda}$ in (4.4). First, the valuation $m/q$ of $y/(\ln x)^{\log_b \lambda}$ is in $q_{\lambda}^{-1}\mathbb{Z}$, and can thus be written $m/q = \tilde{m}/q_{\lambda}$ for $\tilde{m} \in \mathbb{Z}$. Second, $f(x^{1/q}) = y/((\ln x)^{\log_b \lambda} x^{\tilde{m}/q_{\lambda}})$ is in $\mathbb{L}((x^{1/q_{\lambda}}))$, so that
$h(x) := f(x^{q_{\lambda}/q})$ is in $\mathbb{L}((x))$. From $Mf/f = g \in \mathbb{L}(x)$ follows the membership of $g(x^{q_{\lambda}/q}) = Mh/h$ in $\mathbb{L}(x^{q_{\lambda}/q}) \cap \mathbb{L}((x))$. Write $q_{\lambda}/q = n/d$ in lowest terms and observe that $g$ is a series in $x^d$. Define $\tilde{g}(x) := g(x^{n/d})$, a rational series satisfying $\tilde{g}(0) = 1$. Define $\tilde{f}(x) := \prod_{k \ge 0} 1/\tilde{g}(x^{b^k})$ and observe $f(x^{1/q}) = \tilde{f}(x^{1/q_{\lambda}})$. We have obtained a new expression of $y$ of the form (4.4) with $(q, m, g)$ replaced with $(q_{\lambda}, \tilde{m}, \tilde{g})$. In particular, we have made $q$ independent of $y$.

Pick an element $y_0 \in \mathfrak{H}_{\ne 0}$ and obtain its decomposition of the form (4.4), thus fixing $(q, m, g)$ to some $(q_{\lambda}, m_0, g_0)$. Replace $y_0$ with $y_0/c$ so as to set $c$ to 1. Since $\mathfrak{H}_{\ne 0}$ is an $\mathbb{L}(x^{1/*})$-similarity class, the set $V := \{y/y_0 : y \in \mathfrak{H}\}$ is an $\mathbb{L}$-subspace of $\mathbb{L}(x^{1/*})$, and thus, by the previous paragraph, a finite-dimensional $\mathbb{L}$-subspace of $\mathbb{L}(x^{1/q_{\lambda}})$. Let $a(x^{1/q_{\lambda}})$ be the least common denominator of the elements of $V$, with valuation coefficient normalized to one. Write $a(x^{1/q_{\lambda}}) = x^{v/q_{\lambda}} \tilde{a}(x^{1/q_{\lambda}})$ for $\tilde{a}$ of valuation zero. Then $a(x^{1/q_{\lambda}})V$ is a finite-dimensional subspace of $\mathbb{L}[x^{1/q_{\lambda}}]$; as such, it admits a basis $(\tilde{p}_1(x^{1/q_{\lambda}}), \dots, \tilde{p}_s(x^{1/q_{\lambda}}))$ in reduced echelon form w.r.t. ascending degree. The polynomials $\tilde{p}_1, \dots, \tilde{p}_s$ are coprime by the minimality of $a$. Setting $g = g_0a/Ma$ and $p_j = x^{m_0-v}\tilde{p}_j$, the expressions (4.3) form a basis of $\mathfrak{H}$ that satisfies the conditions 2 and 3.

We already have noticed that $\lambda$ is uniquely determined by $\mathfrak{H}$. Suppose that $\mathfrak{H}$ admits a second basis of the form (4.3), with parameters $g', p'_1, \dots, p'_s$ also satisfying 2 and 3. Let $f' = \prod_{k \ge 0} 1/g'(x^{b^k})$. Then $(p_1, \dots, p_s)$ and $(p'_1 f'/f, \dots, p'_s f'/f)$ are two bases of the same $\mathbb{L}$-vector space. In particular, the $p'_j f'/f$ are Laurent polynomials. A Bézout relation $\sum_j u_j p'_j = 1$ for Laurent polynomials $u_j$ then implies that $f'/f = \sum_j u_j (p'_j f'/f)$ is a Laurent polynomial. By symmetry, $f/f'$ is also Laurent, and in fact equal to 1 because $g(0) = g'(0) = 1$. Then the uniqueness of the reduced echelon form of a given vector space implies $(p_1, \dots, p_s) = (p'_1, \dots, p'_s)$. We have thus proved the uniqueness of the tuple $(\lambda, g, p_1, \dots, p_s)$ under conditions 1, 2, and 3. $\square$

A direct consequence of Proposition 4.3 is that for a given $(L)$, $\mathbb{L}(x^{1/*})$-similarity classes are uniquely identified by the pairs $(\lambda, g)$ extracted from the tuples $(\lambda, g, p)$ of the proposition. Indeed, if two classes share the same $(\lambda, g)$, then, in view of the form (4.3) of the bases associated with the two classes, any element of one class must be $\mathbb{L}(x^{1/*})$-similar to any element of the other class. These two classes must therefore be equal. This leads us to the following definition that the reader will compare with Definition 3.21.

**Definition 4.4.** Given a nontrivial $\mathbb{L}(x^{1/*})$-similarity class of $\mathbb{L}(x^{1/*})$-hyper-geometric solutions of $(L)$, let $(\lambda, g, (p_j)_{1 \le j \le s})$ be the uniquely determined tuple of Proposition 4.3. Then the similarity class will be denoted $\mathfrak{H}_{\mathbb{L},\lambda,g}$. We will also write $\mathfrak{R}_{\mathbb{L},\lambda,g}$ for the image of $\mathfrak{H}_{\mathbb{L},\lambda,g \ne 0}$ under $y \mapsto My/y$.

Summing up the previous discussion, the full solution set of the Riccati equation can be described as follows.

**Theorem 4.5.** The solution set $\mathfrak{R}_{\mathbb{L}(x^{1/*})}$ of the Riccati equation (R) in $\mathbb{L}(x^{1/*})$ is a disjoint union

$$ \mathfrak{R}_{\mathbb{L}(x^{1/*})} = \coprod_{(\lambda,g)} \mathfrak{R}_{\mathbb{L},\lambda,g} \qquad (4.5) $$
indexed by a finite set of pairs $(\lambda, g)$, where $\lambda$ is in $\mathbb{L} \cap \Lambda'$ and $g$ is a rational function satisfying $g(0) = 1$, such that the $\mathfrak{H}_{\mathbb{L},\lambda,g}$ partition $\mathfrak{H}_{\mathbb{L},\lambda}$. Each $\mathfrak{R}_{\mathbb{L},\lambda,g}$ can be parametrized bijectively by

$$u(x) = \lambda \frac{\sum_{j=1}^{s} c_j p_j(x^{b/q_{\lambda}})}{\sum_{j=1}^{s} c_j p_j(x^{1/q_{\lambda}})} g(x^{1/q_{\lambda}}) \quad \text{for} \quad (c_1 : \dots : c_s) \in \mathbb{P}^{s-1}(\mathbb{L})$$

where the $p_j$ are those in the uniquely determined tuple of Proposition 4.3. The dimensions $s-1$ of the projective spaces add up to a number that is at most the order $r$ of the linear equation (L) minus the cardinality of the set of pairs $(\lambda, g)$ indexing the disjoint union (4.5).

*Proof.* The classes $\mathfrak{H}_{\lambda,0}$ of Theorem 2.9, applied to $D = \mathbb{L}((x^{1/*}))$ and $F = \mathbb{L}(x^{1/*})$, are the $\mathfrak{H}_{\mathbb{L},\lambda,g}$ of Definition 4.4. The result then follows, using Remark 2.10 and the explicit basis (4.3) of $\mathfrak{H}_{\mathbb{L},\lambda,g}$ provided by Proposition 4.3. $\square$

## 5. Degree bounds for rational solutions

In the present section, we derive degree bounds for rational solutions of the Riccati equation in terms of its order $r$ and coefficient degree $d$. Such bounds will be useful in §8.1 to select relevant Hermite-Padé approximants.

The following easy consequence of Lemma 3.1(c) in (CDDM 2018) will be used throughout; we prove it for self-containedness.

**Lemma 5.1.** Given a field $K$, $P$ and $Q$ in $K[x]$, and any $j \in \mathbb{N}$, $P$ and $Q$ are coprime if and only if so are $M^jP$ and $M^jQ$.

*Proof.* If $P$ and $Q$ are coprime, there exist polynomials $A$ and $B$ such that the relation $AP + BQ = 1$ holds. Therefore, $M^jA M^jP + M^jB M^jQ = 1$, so that $M^jP$ and $M^jQ$ are coprime. Conversely, if $M^jP$ and $M^jQ$ are coprime, then there exist polynomials $A$ and $B$ such that $AM^jP + BM^jQ = 1$. If some $G \in K[x]$ divides both $P$ and $Q$, then $M^jG$ divides 1, so $G$ must be a constant. Therefore, $P$ and $Q$ are coprime. $\square$

**Proposition 5.2.** Assume that $P/Q \in \bar{\mathbb{K}}(x)$, with $\gcd(P,Q) = 1$, is a solution of (R). Then, the following degree bounds hold:

$$\deg P \le B_{\text{num}} := \begin{cases} 2d & (b=2), \\ 4d/b^{r-1} & (b \ge 3), \end{cases} \quad (5.1)$$

$$\deg Q \le B_{\text{den}} := \begin{cases} 2(1 - 1/2^r)d & (b=2), \\ 3d/b^{r-1} & (b \ge 3). \end{cases} \quad (5.2)$$

*Proof.* By Lemma 5.1 over $K = \bar{\mathbb{K}}$, the assumption $\gcd(P,Q) = 1$ implies $\gcd(M^{r-1}P, M^{r-1}Q) = 1$, that is, $M^{r-1}P$ and $M^{r-1}Q$ are coprime. For $n \in \mathbb{N}$,
let $y^{\bar{n}}$ denote the "rising factorial" $\prod_{i=0}^{n-1} M^i y$. In particular, $y^{\bar{0}} = 1$. Multiplying (R) by $Q^{\bar{r}}$ yields

$$ \sum_{i=0}^{r} \ell_i P^{\bar{i}} (M^i Q)^{\bar{r}-i} = 0. \qquad (5.3) $$

Since the factor $M^{r-1}Q$ appears in all terms on the left-hand side of (5.3) but the term for $i = r$, and since $M^{r-1}Q$ and $M^{r-1}P$ are coprime, $M^{r-1}Q$ must divide $\ell_r P^{\bar{r}-1}$. Hence, since $\deg \ell_r \le d$, the degrees of $Q$ and $P$ are related by the inequality

$$ b^{r-1} \deg Q \le d + \frac{b^{r-1}-1}{b-1} \deg P. \qquad (5.4) $$

Furthermore, for (5.3) to hold, at least one term of index $i < r$ must have a degree greater than or equal to the degree of the term of index $r$. Simplifying by the common factor $P^{\bar{i}}$ yields

$$ \deg (\ell_r (M^i P)^{\bar{r}-i}) \le \deg (\ell_i (M^i Q)^{\bar{r}-i}), $$

from which follows

$$ \deg P - \deg Q \le \frac{\deg \ell_i - \deg \ell_r}{(b^r - b^i)/(b-1)} \le \frac{d}{b^{r-1}}. $$

The above inequality yields an upper bound on $\deg P$,

$$ \deg P \le \frac{d}{b^{r-1}} + \deg Q. \qquad (5.5) $$

Substituting it into (5.4), we find

$$ b^{r-1} \deg Q \le d + \frac{b^{r-1}-1}{b-1} \left( \frac{d}{b^{r-1}} + \deg Q \right) = d \frac{b^r - 1}{(b-1)b^{r-1}} + \frac{b^{r-1}-1}{b-1} \deg Q. $$

Multiplying by $b-1$ yields

$$ (b^r - 2b^{r-1} + 1) \deg Q \le \frac{d}{b^{r-1}} (b^r - 1) \le db. \qquad (5.6) $$

For the case $b \ge 3$, using the last term, $db$, in (5.6), and the inequality $b^r \le 3(b^r - 2b^{r-1} + 1)$ yields (5.2). For the case $b = 2$, specializing the left inequality in (5.6) at $b = 2$ yields precisely (5.2). In both cases, (5.1) is a direct consequence of (5.5) and (5.2). $\square$

**Remark 5.3.** The elementary result of Proposition 5.2 is remarkable as it provides a uniform polynomial bound in terms of scalar parameters describing the size of the equation, namely its order $r$ and degree $d$. This is unreachable in the shift operator case, where the degrees of rational solutions of Riccati equations may be exponential in the bit size of the equation. For example, the Charlier polynomials

$$ C_n(x, a) = {}_2F_0 \left( \begin{array}{c} -n, -x \\ - \\ -a^{-1} \end{array} \right) $$
viewed as functions of $x$ with parameters $n \in \mathbb{N}$ and $a > 0$ satisfy the recurrence equation (NIST 2010, §18.22)

$$aC_n(x+1, a) - (a+x)C_n(x, a) + xC_n(x-1, a) + nC_n(x, a) = 0$$

whose degrees of coefficients are bounded by $d = 1$ and whose last coefficient $n$ has bit size $\log n$, while $C_n(x, a)$ has degree $n$. As a result the associated Riccati equation has solutions $C_n(x+1, a)/C_n(x, a)$ whose numerator and denominator cannot be bounded solely as a function of $d$. The fact that the Mahler operator drastically increases the degree is sometimes a difficulty, but here it allows us to obtain the bounds in Lemma 5.2.

*Remark 5.4.* The exponents of $b$ and $d$ in the bounds (5.1) and (5.2) are tight, as we show now. For a given nonzero rational function $u = P/Q$ in its lowest terms with both $P$ and $Q$ monic, let us consider any operator of order $r \ge 2$

$$L = \left( \sum_{k=0}^{r-1} c_k M^k \right) (QM - P) \\ = c_{r-1}(M^{r-1}Q)M^r + \sum_{k=1}^{r-1} (c_{k-1}M^{k-1}Q - c_kM^kP)M^k - c_0P,$$

with constant coefficients $c_0, \dots, c_{r-1}$ in $\overline{\mathbb{K}}$ satisfying $c_{r-1} \ne 0$. The operator admits $M - u$ as a right-hand factor, hence $u$ is a solution of the Riccati equation according to Lemma 2.1. The leading coefficient of $L$ has degree $b^{r-1} \deg Q$. In case $\deg Q \ge b \deg P$, all other coefficients have lower degree. In case $\deg Q < b \deg P$, the coefficient of $M^{r-1}$ has degree $b^{r-1} \deg P$ and the remaining ones have lower degree. Hence, in both cases the maximum degree of the coefficients of $L$ is equal to $d = b^{r-1} \max(\deg P, \deg Q)$. Upon choosing $b \ge 3$ and $\deg P = \deg Q$, the bounds (5.1) and (5.2) become $B_{\text{num}} = 4 \deg P$ and $B_{\text{den}} = 3 \deg Q$, showing that they overshoot by no more than a constant factor 4.

# Part II: Algorithm by exploration of the possible singularities

## 6. Mahlerian variant of Petkovšek's algorithm

In this section, we present an algorithm for computing hypergeometric solutions of linear Mahler equations adapted from Petkovšek's algorithm for difference equations in the usual shift operator. In doing so, we generalize to arbitrary order a previous adaptation of Petkovšek's algorithm to Mahler equations of order 2 due to Roques.

Equations are given with polynomial coefficients over a field $\mathbb{K}$, and we solve them for solutions with coefficients in a field $\mathbb{L}$ satisfying $\overline{\mathbb{K}} \supseteq \mathbb{L} \supseteq \mathbb{K}$. As
a consequence of §3.2, particularly the definition (3.10) for the ramification bound $q_λ$, we have a bound

$$q_{\mathbb{L}} = \operatorname{lcm}_{\lambda \in \mathbb{L} \cap \Lambda} q_{\lambda} \qquad (6.1)$$

on the ramification orders of solutions in $K(x^{1/*})$ of the Riccati equation, so that computing all ramified rational solutions reduces to computing the plain rational solutions of a modified equation. Additionally, the method to be used in §6.2 and §6.3 requires a supplementary ramification in its intermediate calculations: whatever the target ramification order, our method requires the working ramification order to be multiplied by a factor $b^{r-1}$.

### 6.1. Petkovšek's classical algorithm

In the case of the linear classical difference equation

$$\ell_r(x)y(x+r) + \cdots + \ell_0(x)y(x) = \sum_{i=0}^{r} \ell_i(x)y(x+i) = 0 \qquad (6.2)$$

and the corresponding Riccati equation

$$\ell_r(x)u(x) \cdots u(x+r-1) + \cdots + \ell_1(x)u(x) + \ell_0(x) = \sum_{i=0}^{r} \ell_i(x) \prod_{j=0}^{i-1} u(x+j) = 0, \qquad (6.3)$$

an algorithm due to Petkovšek (1992) is known to solve (6.3) for all its rational function solutions, or equivalently to find all first-order right-hand factors of (6.2).

The algorithm is based on the concept of a *Gosper-Petkovšek form*: for any rational function $u(x) \in \mathbb{L}(x)_{\neq 0}$, there exist a constant $\zeta \in \mathbb{L}_{\neq 0}$ and monic polynomials $A(x), B(x), C(x)$ in $\mathbb{L}[x]$ satisfying:

1. $u(x) = \zeta \frac{C(x+1)}{C(x)} \frac{A(x)}{B(x)}$,

2. $A(x)$ and $C(x)$ are coprime,

3. $B(x)$ and $C(x+1)$ are coprime,

4. $A(x)$ and $B(x+i)$ are coprime for all $i \ge 0$.

Now, any potential nonzero rational solution $u$ of the Riccati equation (6.3) leads to the necessary relation

$$\sum_{i=0}^{r} \ell_i(x) \zeta^i C(x+i) \left( \prod_{j=0}^{i-1} A(x+j) \right) \left( \prod_{j=i}^{r-1} B(x+j) \right) = 0. \qquad (6.4)$$

Here, $A(x)$ appears in all the terms of the sum but the one for $i=0$. As it is coprime to all forward shifts of $B(x)$ and to $C(x)$, it must divide $\ell_0(x)$. Similarly, $B(x+r-1)$ appears in every term but the one for $i=r$, and so must divide $\ell_r(x)$. This motivates iterating on all pairs of monic divisors of $\ell_0$ and $\ell_r$. Given monic $A(x) | \ell_0(x)$ and $B(x) | \ell_r(x-r+1)$, where divisibility is meant in $\mathbb{L}[x]$, the
leading coefficient of the left-hand side of (6.4) is independent of the choice of a monic $C(x)$. So this leading coefficient yields an algebraic equation in $\zeta$. Each choice of a solution $\zeta \in \mathbb{L}$ turns (6.4) into an equation that can be solved for polynomial solutions $C(x) \in \mathbb{L}[x]$. The algorithm then gathers and returns all found $(\zeta, A(x), B(x), C(x))$.

Remark that the classical definition of the literature, quantified over all integers $i \in \mathbb{N}$ and reproduced as point 4 above, is stronger than needed: the proof has used the property for $0 \le i \le r - 1$ only. In the Mahler case, it will prove important to use the weaker constraint to get an algorithm. The reader should also compare point 4 of the present section with its analogues in §6.2 and §6.3.

## 6.2. Roques's algorithm for order 2

Inspired by Hendriks' works for usual difference equations (Hendriks 1998) and for $q$-difference equations (Hendriks 1997), Roques (2018, §6.2) recently presented an analogue of Petkovšek's algorithm for Mahler equations of order 2, that is, equations of the form

$$ \ell_2(x)u(x)u(x^b) + \ell_1(x)u(x) + \ell_0(x) = 0. \quad (6.5) $$

Roques's original presentation determines along his algorithm a suitable extension of $\mathbb{K}$ sufficient to obtain all solutions from $\bar{\mathbb{K}}(x^{1/*})$. Here, we present a variant that computes with an input field $\mathbb{L}$, so as to be consistent with the rest of our text.

After finding a bound $q$ such that all solutions $u \in \mathbb{L}(x^{1/*})$ are in fact in $\mathbb{L}(x^{1/q})$, let us introduce a new indeterminate $t$ for which $x = t^{qb}$, so as to prove, for any given $u \in \mathbb{L}(x^{1/*})$, the existence of a nonzero constant $\zeta \in \mathbb{L}_{\ne 0}$ and monic polynomials $A, B, C \in \mathbb{L}[t]$ satisfying

$$ 1. u(x) = \zeta \frac{C(x^{1/q}) A(x^{1/qb})}{C(x^{1/qb}) B(x^{1/qb})}, \text{ that is, } u(t^{qb}) = \zeta \frac{C(t^b) A(t)}{C(t) B(t)}, $$

2. $A$ and $C$ are coprime,

3. $B$ and $MC$ are coprime,

4. $A$ and $M^iB$ are coprime for all $i \in \{0, 1\}$,

5. $C$ and $MC$ are coprime,

where $M$ acts on $\mathbb{L}[t]$ by substituting $t^b$ for $t$.

Now, any nonzero ramified rational solution $u$ of the Riccati equation (6.5) leads to the necessary relation

$$ \ell_2(t^{qb})\zeta^2 C(t^{b^2})A(t)A(t^b) + \ell_1(t^{qb})\zeta C(t^b)A(t)B(t^b) + \ell_0(t^{qb})C(t)B(t)B(t^b) = 0. \quad (6.6) $$

This time, one finds that $A(t)$ must divide $\ell_0(t^{qb})$, while $B(t^b)$ must divide $\ell_2(t^{qb})$, implying that $B(t)$ must divide $\ell_2(t^q)$.
Roques's approach continues in a way similar in spirit to Petkovšek's, although technical reasons require to exchange the order of steps. Roques indeed finds a linear constraint on the degree of $C$ by putting apart the term with $\zeta^2$ in (6.6). Finding the coefficients of $C$ then amounts to solving a finite linear system. As a side remark, in Roques's text valuations and degrees are treated interchangeably, and the calculation of the coefficient $\zeta$ is obscured. We will clarify this in the next section.

### 6.3. A new algorithm for higher-order Mahler equations

In this section, we deal with Riccati Mahler equations (R) of general order $r \ge 2$. The core of the section describes an algorithm for solving for (plain) rational functions from $\mathbb{L}(x)$.

In order to go beyond order 2, we observe that point 4 in Roques's definition requires a coprimality only for $0 \le i \le r-1 = 1$ (when $r=2$), whereas point 4 in the shift situation requires a coprimality for $0 \le i$, although the proof in §6.1 for an equation of order $r$ only uses the cases $0 \le i \le r-1$. This motivates us to introduce the following notion of a bounded Gosper-Petkovšek form.

**Definition 6.1.** Given a rational function $u(x) \in \mathbb{L}(x)_{\ne 0}$ and an integer $r \ge 2$, a bounded Gosper-Petkovšek form of order $r$ for $u(x)$ is a tuple $(\zeta, A, B, C) \in \mathbb{L}_{\ne 0} \times \mathbb{L}[t]^3$, with $A, B, C$ monic polynomials in a new indeterminate $t$, such that:

$$1. \ u(t^{b^{r-1}}) = \zeta \frac{C(t^b)}{C(t)} \frac{A(t^{b^{r-1}})}{B(t)},$$

$$2. \ M^{r-1}A \text{ and } C \text{ are coprime},$$

$$3. \ B \text{ and } MC \text{ are coprime},$$

$$4. \ M^i A \text{ and } B \text{ are coprime for all } i \in \{0, \dots, r-1\},$$

where $M$ acts on $\mathbb{L}[t]$ by substituting $t^b$ for $t$.

By Lemma 5.1 over $K = \mathbb{L}$, point 4 above can be restated equivalently into: $M^{r-1}A$ and $M^i B$ are coprime for all $i \in \{0, \dots, r-1\}$. Also remark that for a polynomial triple $(A, B, C)$ satisfying our definition for $r=2$, the polynomial triple $(MA, B, C)$ satisfies Roques's definition (with $q=1$), at least provided $\gcd(C, MC) = 1$. A constructive proof of the existence of bounded Gosper-Petkovšek forms will be provided in §6.4.

Now, consider any potential nonzero rational solution $u$ of the Riccati equation (R), represented by one of its bounded Gosper-Petkovšek forms. Substituting the Gosper-Petkovšek form for $u$ and $t^{b^{r-1}}$ for $x$ in (R) leads after canceling denominators to the necessary relation

$$\tilde{L}(t, \zeta M) C = 0, \qquad (6.7)$$

where $\tilde{L}$ is the operator in $\mathbb{L}[t]\langle M \rangle$ defined by

$$\tilde{L}(t, M) = \sum_{k=0}^{r} \left( M^{r-1} \ell_k(t) \times \prod_{j=0}^{k-1} M^{r-1+j} A(t) \times \prod_{j=k}^{r-1} M^j B(t) \right) M^k. \quad (6.8)$$
In this formula, note that beside the polynomials A and B that are by definition in the ring $\mathbb{L}[t]$, we have noted $\ell_k(t)$ for the result of the substitution of $t$ for $x$ in $\ell_k$. The factor $M^{r-1}A$ appears in all the terms of the expansion of (6.7), except for the one corresponding to $k=0$,

$$M^{r-1}\ell_0(t) \times C(t) \times \prod_{j=0}^{r-1} M^j B(t);$$

using points 2 and 4 in Definition 6.1, $M^{r-1}A$ must divide $M^{r-1}\ell_0$, and Lemma 5.1 implies that $A$ must divide $\ell_0$. Similarly, $M^{r-1}B$ appears in all the terms of (6.7), except for the one corresponding to $k=r$,

$$M^{r-1}\ell_r(t) \times \zeta^r M^r C(t) \times \prod_{j=0}^{r-1} M^{r-1+j} A(t);$$

thus, $M^{r-1}B$ must divide $M^{r-1}\ell_r$, and Lemma 5.1 implies that $B$ must divide $\ell_r$.

To complete the pairs $(A, B)$ into candidate tuples $(\zeta, A, B, C)$ delivering rational functions in $\mathbb{L}(t)$, we interpret (6.7) as an analogue of (3.9) in Proposition 3.8. To this end, we start by introducing in analogy with the construction of $\mathfrak{D}$ in §3.1 a ring whose coefficients are series in $t^{-1}$, in the form

$$\mathfrak{D}' = \bigoplus_{\lambda \in \bar{\mathbb{K}} \setminus \mathbb{K}_0} (\ln(t^{-1}))^{\log_b \lambda} \bar{\mathbb{K}}((t^{-1})^{1/*})).$$

Here, the $(\ln(t^{-1}))^{\log_b \lambda}$ denote elements of some new family $(e'_\lambda)_{\lambda \in \bar{\mathbb{K}} \setminus \mathbb{K}_0}$) of generators satisfying the algebraic relations $e'_\lambda e'_{\lambda'} = e'_{\lambda\lambda'}$ and $e'_1 = 1$. The theory of $\mathfrak{D}'$ parallels that of $\mathfrak{D}$, including the results developed in §3.2 to find the possible logarithmic parts $(\ln x)^{\log_b \lambda}$ and valuations of the corresponding Puiseux series in a solution of (L) by using the lower Newton polygon.

**Definition 6.2.** Let $L$ be a general linear Mahler operator. The upper convex hull of the set of points $(b^k, j) \in \mathbb{R}^2$ of Definition 3.4 is called the upper Newton polygon (of $L$). In analogy with Definition 3.6, we index its edges by decreasing slopes and denote by $\xi_j(X)$ the characteristic polynomial of its $j$th edge. Define $Z(L)$ to be the union of the sets of roots in $\bar{\mathbb{K}}$ of the $\xi_j$,

$$Z(L) = \{\zeta \in \bar{\mathbb{K}} : \exists j, \xi_j(\zeta) = 0\}. \qquad (6.9)$$

The $j$th edge will be called $\zeta$-admissible (for $L$) if $\xi_j(\zeta) = 0$.

Now, an analogue of Proposition 3.8 holds, with the difference that the parameter $\lambda$ of the new proposition needs to be the root of some characteristic polynomial associated with an edge of an upper Newton polygon, rather than a lower Newton polygon. Consequently, the constant $\zeta$ in (6.7) must be an element of $Z(\bar{L})$. In this situation, the degree of $C$ in $t$ must be the opposite of the slope of some edge associated with the upper Newton polygon: the degree of a
solution $C$ of (6.7) in $t$ is the opposite of the valuation with respect to $t^{-1}$ of $C$, which is given by the lower Newton polygon of $\tilde{L}(t^{-1}, M)$; the upper Newton polygon of $\tilde{L}(t, M)$ is the symmetric with respect to a horizontal axis of the lower Newton polygon of $\tilde{L}(t^{-1}, M)$; edges associated by this symmetry have opposite slopes; the result follows.

We summarize the previous discussion in the following result.

**Proposition 6.3.** Let $(\zeta, A, B, C)$ be a bounded Gosper-Petkovšek form of order $r$ of any rational solution $u$ of (R), and let $\tilde{L}$ be defined by (6.8). Then, necessarily:

* $\zeta \in Z(\tilde{L})$,

* $A$ divides $\ell_0(t)$,

* $B$ divides $\ell_r(t)$,

* $C$ is a nonzero polynomial solution of (6.7), implying $\deg C$ is the opposite of a nonpositive integer slope of an edge associated with the upper Newton polygon of $\tilde{L}$.

The previous considerations lead to Algorithm 1, whose general structure is the following:

* a loop over candidates $(A, B)$ is set up at step (B) from the lower Newton polygon of the input $L$;

* candidates $\zeta$, then degrees for candidates $C$ at step (B)(3)(a), are obtained from the upper Newton polygon of the auxiliary operator $\tilde{L}$;

* solving for $C$ at step (B)(3)(b)(i) is done by appealing to our algorithm for polynomial solutions of bounded degree in (CDDM 2018, §2.6, Algorithm 5);

* to avoid redundancy in the output, the cleaning step (C) makes sure to return a partition by enforcing that no parametrization is included in another.

Regarding the last item, the cleaning step (C) removes redundant elements of $\mathcal{W}$ by a simple loop: while $\mathcal{W}$ contains two distinct elements $u^{(1)}$ and $u^{(2)}$ with inclusion of $u^{(1)}$ in $u^{(2)}$, it removes $u^{(1)}$ from $\mathcal{W}$. This is effective provided we have an algorithmic inclusion test available. To describe one, observe that, given two parametrized rational functions, $u^{(1)}$ with parameters $c^{(1)} = (c_1^{(1)}, \dots, c_{s_1^{(1)}})$ and $u^{(2)}$ with parameters $c^{(2)} = (c_1^{(2)}, \dots, c_{s_2^{(2)}})$, inclusion of $u^{(1)}$ in $u^{(2)}$ is only possible if $s_1 \le s_2$, in which case getting rid of denominators in the equation $u^{(1)} = u^{(2)}$ and equating like powers of $x$ results in a linear system in $c^{(2)}$ linearly parametrized by $c^{(1)}$. Either this system has no nonzero solution, proving noninclusion, or it provides a parametrization proving inclusion. We thus have an algorithm by linear algebra over $\mathbb{L}$.

We now prove a first part of the correctness of the algorithm: no rational solution is lost.
Input: A Riccati Mahler equation (R) with coefficients $\ell_k(x) \in \mathbb{K}[x]$. Some intermediate field $\mathbb{L}$, that is, a field satisfying $\bar{\mathbb{K}} \supseteq \mathbb{L} \supseteq \mathbb{K}$.

Output: The set of rational functions $u \in \mathbb{L}(x)$ that solve (R).

(A) Set $\mathcal{W} := \emptyset$.

(B) For each monic factor $A(t) \in \mathbb{L}[t]$ of $\ell_0(t)$, for each monic factor $B(t) \in \mathbb{L}[t]$ of $\ell_r(t)$ such that $M^i A$ and $B$ are coprime for $0 \le i < r$:

(1) compute $\tilde{L}(t, M)$ by (6.8),

(2) compute the upper Newton polygon of $\tilde{L}$, the set $Z(\tilde{L}) \cap \mathbb{L}$ of roots $\zeta$ in $\mathbb{L}$ of the associated characteristic polynomials,

(3) for each $\zeta$ in $Z(\tilde{L}) \cap \mathbb{L}$:

(a) compute the maximum $\Delta_\zeta$ of the integer values of the opposites of the slopes of the $\zeta$-admissible edges (for $\tilde{L}$),

(b) if $\Delta_\zeta \ge 0$:

(i) compute a basis $(C_i)_{1 \le i \le s}$ of solutions in $\mathbb{L}[t]_{\le \Delta_\zeta}$ of the equation $\tilde{L}(t, \zeta M) C = 0$,

(ii) if $s > 0$:

(α) set $C := \sum_{i=1}^s c_i C_i$ for formal parameters $c_i$,

(β) normalize the rational function

$$ \tilde{u}(t) := \zeta \frac{C(t^b) A(t^{b^{r-1}})}{C(t) B(t)}, $$

which is an element of $\mathbb{L}(c_1, \dots, c_s)(t^{b^{r-1}})$, so as to identify $u(x) \in \mathbb{L}(c_1, \dots, c_s)(x)$ such that $u(t^{b^{r-1}}) = \tilde{u}(t)$,

(γ) augment $\mathcal{W}$ with $u$.

(C) Remove redundant elements from $\mathcal{W}$ by the method described before Proposition 6.4 and return the resulting $\mathcal{W}$.

**Algorithm 1:** Rational solutions to a Riccati Mahler equation. Compare with the efficiency improvements in Algorithm 3.

**Proposition 6.4.** *Algorithm 1 computes the set of rational solutions of the Riccati equation (R) as a union of sets parametrized by finite-dimensional $\mathbb{L}$-projective spaces.*

*Proof.* Observe that Algorithm 1 iterates on all tuples $(\zeta, A, B, \Gamma(c))$ in $\mathbb{L}_{\neq 0} \times \mathbb{L}[t]^2 \times \mathbb{L}(c_1, c_2, \dots)[t]$ such that, for any values of the parameters $c_i$ in $\mathbb{L}$, the tuple $(\zeta, A, B, C)$ where $C = \Gamma(c)$ satisfies: the necessary conditions (i), (ii), and (iii) of Proposition 6.3, point 4 of Definition 6.1, a degree bound on $C$ that is necessarily satisfied by the polynomial kernel of $\tilde{L}(t, \zeta M)$, cf. the necessary condition (iv) in the same proposition. The algorithm therefore represents a set that is less constrained that the set of bounded Gosper-Petkovšek forms of order $r$ of $u$ for any rational solution $u$ of (R); in particular, it represents them all, and, by the existence of bounded Gosper-Petkovšek forms (see §6.4), the algorithm must find any solution $u$ of (R). Conversely, for any $u(x)$ element of the output $\mathcal{W}$, obtained by the algorithm from a tuple $(\zeta, A, B, C)$, expanding (6.7)
by using (6.8), then dividing by the relevant product shows that $u(x)$ solves (R). We have thus proved that Algorithm 1 computes all rational solutions. $\square$

For the end of the present section, we call a block the image of a similarity class of nonzero hypergeometric elements under $y \mapsto My/y$. The following theorem states the correctness of Algorithm 1.

**Theorem 6.5.** Algorithm 1 computes the set of rational solutions of the Riccati equation (R) as a disjoint union of sets bijectively parameterized by finite-dimensional $\mathbb{L}$-projective spaces. Equivalently, the output $\mathcal{U}$ of Algorithm 1 consists of formal parametrizations of blocks, with exactly one parametrization for each block, and each such parametrization is a bijection.

*Proof.* Each element $U(c) = (c_1, \dots, c_s)$ of $\mathcal{U}$ can be interpreted as a parametrization

$$ \gamma = (\gamma_1 : \dots : \gamma_s) \in \mathbb{P}^{s-1}(\mathbb{L}) \mapsto U(\gamma) = \zeta \frac{MC M^{r-1} A}{C B} \in \mathbb{L}(x). \quad (6.10) $$

Fix $y_0$, a nonzero solution of $My_0 = \zeta \frac{M^{r-1}A}{B}y_0$. The solutions of $My = U(\gamma)y$ when $\gamma$ ranges in $(\mathbb{L}^s)_{\neq 0}$ form a vector space $\mathfrak{G}$ that is exactly $C(\mathbb{L}^s)y_0$. Because $C$ is a polynomial in $x$, the set $\mathfrak{G}_{\neq 0}$ is included in a similarity class.

For each similarity class $\mathfrak{H}_{\neq 0}$, the finite-dimensional space $\mathfrak{H}$ is by Proposition 6.4 covered as a finite union of vector spaces $\mathfrak{G}$ obtained from elements $U \in \mathcal{U}$. Let $d$ denote the dimension of $\mathfrak{H}$. Assume that no $\mathfrak{G}$ has dimension $d$. Then, each $\mathfrak{G}$ is in some hyperplane defined by some nonzero linear form $\phi$. Fix a coordinate system of $\mathfrak{H}$ and let $S$ be the parametrized curve $t \mapsto (1, t, \dots, t^{d-1})$, whose image is included in $\mathfrak{H}$. The curve can only meet a given $\mathfrak{G}$ at finitely many intersections, provided by the zeros of the polynomial $\phi(S(t))$. Because the curve has infinitely many points, this contradicts that $\mathfrak{H}$ is covered. So at least one of the $\mathfrak{G}$ is equal to $\mathfrak{H}$.

By the absence of redundancy enforced at step (C) in the algorithm, no two distinct elements of $\mathcal{U}$ can produce vector spaces $\mathfrak{G}_1$ and $\mathfrak{G}_2$ with $\mathfrak{G}_1 \subseteq \mathfrak{G}_2$. So any vector space $\mathfrak{H}$ of a similarity class $\mathfrak{H}_{\neq 0}$ is obtained exactly once as a $\mathfrak{G}$, and only such $\mathfrak{H}$ are obtained. The result for blocks then follows.

Because the algorithm constructs $C$ in (6.10) as a formal sum $C = \sum_{j=1}^s c_j C_j$ for $\mathbb{L}$-linearly independent polynomials $C_j \in \mathbb{L}[x]$, the $C_j y_0$ form an $\mathbb{L}$-basis of $C(\mathbb{L}^s)y_0$. The implied bijection from $\mathbb{L}^s$ to the suitable $\mathfrak{H}$ induces a bijection from $\mathbb{P}^{s-1}(\mathbb{L})$ to the block associated with $\mathfrak{H}_{\neq 0}$. $\square$

**Corollary 6.6.** The set of $\mathbb{L}(x)$-hypergeometric solutions of the Mahler equation (L) can be described as a union of $\mathbb{L}$-vector spaces $\mathfrak{H}$ in direct sum parametrized by the output $\mathcal{U}$ as follows: for given $U \in \mathcal{U}$, the vector space $\mathfrak{H}$ is spanned by a basis $(y_1, \dots, y_s)$ where each $y_i$ is any nonzero solution of the equation $My_i = U(e_i)y_i$ for the $i$th element of the canonical basis.

*Proof.* The statement is a direct consequence of the last paragraph in the proof of Theorem 6.5. $\square$
*Remark 6.7* (efficiency of Algorithm 1). The double loop over *A* and *B* induces an exponential behavior, making this basic algorithm inefficient. We will discuss several pruning strategies in §7. Furthermore, the algorithm considers tuples (ζ, *A*, *B*, *C*) that provide solutions but are not bounded Gosper–Petkovšek forms. These tuples are redundant because the algorithm also produces the same solutions in bounded Gosper–Petkovšek form. From the point of view of efficiency, the problem is not in rejecting tuples that are not bounded Gosper–Petkovšek forms, but in not computing too many (redundant) candidates in the first place. The necessity to determine the polynomials *C* late in the algorithm makes it impossible to consider only bounded Gosper–Petkovšek forms, and is the main cause for getting repeated solutions in the course of the algorithm (see also the Remark 6.10 below).

*Remark 6.8* (solving large linear systems). When searching for polynomial solutions of the auxiliary equations $\tilde{L}(t, \zeta M)C = 0$, it is crucial for performance to use the fast algorithms introduced in (CDDM 2018). We take the opportunity to clarify a minor point of confusion in that work. In (CDDM 2018, Algorithm 2), we appeal to the linear solving algorithm of Ibarra, Moran, and Hui (1982) for computing the kernel of an $m \times n$ matrix $A$ in $O(m^{\omega-1}n)$ field operations. This algorithm is based on computing an LSP decomposition of $A$, that is, a decomposition as a product of a lower square matrix $L$, a semi-upper triangular matrix $S$, and a permutation matrix $P$. Strictly speaking, however, this approach is only valid when $m \le n$ (fewer rows than columns). As our need is for an $m \times n$ matrix $A$ with $m \ge n$ (more rows than columns), a workaround is to apply the LSP decomposition algorithm to the transpose $A^T$, thus obtaining a PSU decomposition of $A$, that is, as a permutation matrix $P$, a semi-lower triangular matrix $S$, and an upper square matrix $U$. To solve for the (right) kernel, we observe $\ker A = U^{-1}(\ker S)$. As $S$ reduces to a lower triangular matrix with nonzero diagonal elements when the zero columns are deleted, finding its kernel is immediate. Inverting $U$ takes $O(n^\omega)$ operations. So $\ker A$ is obtained in $O(m/n)$ square matrix products, hence in complexity $O(mn^{\omega-1})$. Taking the rank $r$ of $A$ into account and refining the analysis yields the complexity $O(mnr^{\omega-2})$ announced in (CDDM 2018, Prop. 2.14).

To compute all *ramified* rational solutions of (R), we now propose Algorithm 2, which is a simple variant of Algorithm 1. By the property of $q_L$ to be a uniform bound on the ramification order of ramified rational solutions, Algorithm 2 is also correct, in the sense that it satisfies an analogue of Theorem 6.5, where rational solutions are replaced by ramified rational solutions. An obvious adaptation of Corollary 6.6 to $\mathbb{L}(x^{1/*})$-hypergeometric solutions also holds.

## 6.4. Existence and computation of bounded Gosper–Petkovšek forms

The following lemma proves the existence of bounded Gosper–Petkovšek forms at all orders, and its proof provides an algorithm for putting a rational function in bounded Gosper–Petkovšek form of a given order.
Input: A Riccati Mahler equation (R) with coefficients $\ell_k(x) \in \mathbb{K}[x]$. Some intermediate field $\mathbb{L}$, that is, a field satisfying $\mathbb{K} \supseteq \mathbb{L} \supseteq \mathbb{K}$.

Output: The set of ramified rational functions $u \in \mathbb{L}(x^{1/*})$ that solve (R).

(A) Compute the lower Newton polygon of $L := \sum_{k=0}^r \ell_k(x)M^k$, the associated characteristic polynomials, the set $\Lambda \cap \mathbb{L}$, and the ramification bound $q_{\mathbb{L}}$ by (6.1).

(B) Call Algorithm 1 for computing the rational solutions of $L(x^{q_{\mathbb{L}}}, M)$ in $\mathbb{L}(x)$.

(C) Substitute $x^{1/q_{\mathbb{L}}}$ for $x$ in the obtained solutions and return the resulting set.

**Algorithm 2:** Ramified rational solutions to a Riccati Mahler equation.

**Lemma 6.9.** Given any two coprime monic polynomials $P$ and $Q$ in $\mathbb{K}[x]$, define a sequence of triples of polynomials given for all $k \in \mathbb{N}_{\neq 0}$ by

$$ (A_1, B_1, C_1) = (P, Q, 1) \qquad (6.11) $$

$$ (A_{k+1}, B_{k+1}, C_{k+1}) = \left( \frac{A_k}{G_k}, \frac{MB_k}{G_k}, MC_k \times (M^0 G_k \cdots M^{k-1} G_k) \right), \quad (6.12) $$

where $G_k = \gcd(A_k, MB_k)$. Then, for all $k \in \mathbb{N}_{\neq 0}$, $(1, A_k, B_k, C_k)$ is a bounded Gosper-Petkovšek form of order $k$ for the rational function $P/Q$. Additionally, there exists $k \in \mathbb{N}_{\neq 0}$ satisfying

$$ A_{k+i} = A_k, \quad B_{k+i} = M^i B_k, \quad C_{k+i} = M^i C_k, \quad \text{for all } i \in \mathbb{N}. \qquad (6.13) $$

*Proof.* The proof is by induction on $k$. The case $k=1$ is immediate by (6.11). For some $k \in \mathbb{N}_{\neq 0}$, assume that $(1, A_k, B_k, C_k)$ is a bounded Gosper-Petkovšek form of order $k$ for the rational function $P/Q$ and that (6.12) holds. This provides the formulas:

$$ M^{k-1} \frac{P}{Q} = \frac{MC_k M^{k-1} A_k}{C_k B_k}, \qquad (6.14) $$

$$ \gcd(M^{k-1} A_k, C_k) = 1, \qquad (6.15) $$

$$ \gcd(B_k, MC_k) = 1, \qquad (6.16) $$

$$ \gcd(M^i A_k, B_k) = 1 \quad \text{for all} \quad 0 \le i < k. \qquad (6.17) $$

Then, first observe

$$ \frac{MC_{k+1}}{C_{k+1}} = \frac{M^2 C_k}{MC_k} \frac{M^k G_k}{G_k} $$

and apply $M$ to (6.14), so that

$$ M^k \frac{P}{Q} = \frac{M^2 C_k}{MC_k} \frac{M^k (A_{k+1} G_k)}{B_{k+1} G_k} = \frac{MC_{k+1}}{C_{k+1}} \frac{M^k A_{k+1}}{B_{k+1}}. $$

Applying Lemma 5.1 with $K = \mathbb{K}$ and $j = 1$ to (6.15–6.17) implies

$$ \gcd(M^{i+1}(A_{k+1}G_k), B_{k+1}G_k) = 1 \quad \text{for all} \quad 0 \le i < k, \\ \gcd(M^k(A_{k+1}G_k), MC_k) = \gcd(B_{k+1}G_k, M^2C_k) = 1, $$
so that in particular

$$ \gcd(M^{i+1}A_{k+1}, B_{k+1}) = \gcd(M^{i+1}A_{k+1}, G_k) = \gcd(M^{i+1}G_k, B_{k+1}) = 1, \qquad (6.18) $$

$$ \gcd(M^k A_{k+1}, MC_k) = \gcd(B_{k+1}, M^2 C_k) = 1. \qquad (6.19) $$

The construction (6.12) ensures $\gcd(M^0 A_{k+1}, B_{k+1}) = 1$, so combining with the left $gcd$ in (6.18) shows (6.17) at $k+1$. Applying Lemma 5.1 with $j = k - i - 1$ to the middle $gcd$ in (6.18) next yields

$$ \gcd(M^k A_{k+1}, M^{k-i-1} G_k) = 1. $$

Considering those $gcds$ for $0 \le i < k$ as well as the left $gcd$ in (6.19) proves (6.15) at $k+1$. Combining the right $gcd$ in (6.18) for $0 \le i < k$ together with the right $gcd$ in (6.19) proves (6.16) at $k+1$. We have obtained that $(1, A_{k+1}, B_{k+1}, C_{k+1})$ is a bounded Gosper-Petkovšek form of order $k+1$ for the rational function $P/Q$.

Last, in view of (6.12), the existence of $k$ satisfying (6.13) is equivalent to the existence of $k$ such that $G_{k+i} = 1$ for all $i \in \mathbb{N}$. This $k$ exists because $\deg A_k$ cannot decrease indefinitely. $\square$

*Remark 6.10.* Bounded Gosper-Petkovšek forms are not unique, even for a given order. An example with $b = 2$ is given for order $r = 3$ by $(\zeta, A, B, C) = (1, x+1, 1, 1)$ and $(\zeta', A', B', C') = (1, 1, 1, x^4 - 1)$, where the first form is computed by the direct use of the recurrence in Lemma 6.9.

*Remark 6.11.* Point 5 in our presentation of Roques's definition (see §6.2) is not enforced in our definition of bounded Gosper-Petkovšek forms. (Also compare with the constraints on $r$, $s$, and $t$ in (Roques 2018, §6.2, top of page 345).) The mere construction of $C_{k+1}$ in (6.12) makes the coprimality between $C$ and $MC$ generally impossible, at least for $r \ge 3$.

# 7. A more practical algorithm

As experiments show (see §10.3 and especially the timings in Table 2), a direct implementation of Algorithm 1 is inefficient because of the large number of pairs of divisors $(A, B)$ it has to consider. In the present section, we describe several ways to mitigate this issue. The results are rather technical. They are not used in the rest of the article, apart from the fact that the resulting improved variant algorithm is tested beside our other algorithms in §10.

## 7.1. Ensuring coprimality

Given a polynomial $f \in \mathbb{L}[t]$, let $\text{irred}(f)$ denote the set of its monic irreducible factors and $\text{val}_p f$ denote the multiplicity of a monic irreducible polynomial $p$ in $f$. By representing a monic divisor of $\ell_0$ by the family $(\alpha_p)_{p \in \text{irred}(\ell_0)}$ of exponents in its factorization $\prod_{p \in \text{irred}(\ell_0)} p^{\alpha_p}$ into monic irreducible factors, the set of such divisors can be viewed as the Cartesian product

$$ \mathcal{A} = \prod_{p \in \text{irred}(\ell_0)} [0, \text{val}_p \ell_0], $$
where $[a, b]$ denotes the integer interval $\{a, \dots, b\}$. In the same way, the set of monic divisors of $\ell_r$ is represented by

$$ \mathcal{B} = \prod_{q \in \text{irred}(\ell_r)} [0, \text{val}_q \ell_r]. $$

So, the main loop (B) in Algorithm 1 can be viewed as parametrized by a pair $(\alpha, \beta)$ in the product $\mathcal{A} \times \mathcal{B}$. However, we are only interested in pairs $(A, B)$ that satisfy $\gcd(M^s A, B) = 1$ for all $0 \le s < r$, and the basic algorithm has to test all those gcds in the double loop over $A$ and $B$.

To explain how to prune $\mathcal{A} \times \mathcal{B}$, let us imagine that $\ell_r$ has a monic irreducible factor $q$ with multiplicity $b > 0$ and that $\ell_0$ has a monic irreducible factor $p$ with multiplicity $a > 0$ such that $q$ is a factor of some $M^s p$, $0 \le s < r$. In such a situation, we say that the factor $p$ is forbidden (in $A$) by the factor $q$ (of $B$). When a pair $(\alpha, \beta)$ of tuples is in $\mathcal{A} \times \mathcal{B}$, the pair $(\alpha_p, \beta_q)$ of integers lies in the Cartesian product $[0, a] \times [0, b]$. Nonetheless, only the parts $[0, a] \times \{0\}$ and $\{0\} \times [1, b]$ have to be considered, as $(\alpha_p, \beta_q)$ outside of these parts violate the coprimality condition. The previous considerations are independent of the other coordinates $\alpha_{p'}$ and $\beta_{q'}$, so only a fraction $(a+b+1)/((a+1)(b+1))$ of the whole product $\mathcal{A} \times \mathcal{B}$ is useful. Taking several pairs $(p, q)$ into account, we expect only a small fraction of $\mathcal{A} \times \mathcal{B}$ to remain.

We now make this observation algorithmic. Given a monic irreducible factor $q$ of $\ell_r$, we denote by $\mathcal{F}(q)$ the set of monic irreducible factors $p$ of $\ell_0$ forbidden by $q$. Let $\mathcal{R}$ be the subset of irred($\ell_r$) consisting of those factors $q$ of $\ell_r$ such that $\mathcal{F}(q)$ is nonempty. Given a $B$, a monic irreducible factor $q$ of $B$ restricts the choice of the useful $A$ if it is in $\mathcal{R}$, and places no restriction on $A$ if it is not in $\mathcal{R}$. Next, for a subset $\pi \subseteq \mathcal{R}$, let $\mathcal{F}(\pi)$ denote $\bigcup_{q \in \pi} \mathcal{F}(q)$. The set of useful pairs $(A, B)$ can now be seen to be the disjoint union over the subsets $\pi \subseteq \mathcal{R}$ of the Cartesian products

$$ \mathcal{C}_{\pi} = \prod_{p \in \text{irred}(\ell_0)} \mathcal{A}(p, \pi) \times \prod_{q \in \text{irred}(\ell_r)} \mathcal{B}(\pi, q), $$

where

$$ \mathcal{A}(p, \pi) = \begin{cases} [0, \text{val}_p \ell_0] & \text{if } p \notin \mathcal{F}(\pi), \\ \{0\} & \text{if } p \in \mathcal{F}(\pi), \end{cases} \quad \mathcal{B}(\pi, q) = \begin{cases} [0, \text{val}_q \ell_r] & \text{if } q \notin \mathcal{R}, \\ [1, \text{val}_q \ell_r] & \text{if } q \in \pi, \\ \{0\} & \text{if } q \in \mathcal{R} \setminus \pi. \end{cases} $$

Rather than precomputing the $\mathcal{A}(p, \pi)$ and $\mathcal{B}(\pi, q)$ explicitly, in our algorithm we discard useless pairs on the fly.

In the following, we call *factored representation* of a polynomial its representation as a power product of monic irreducible factors. An optimization of Algorithm 1 is obtained by just changing the iteration of loop (B) into

(B) For all $B := \prod_{q \in \text{irred}(\ell_r)} q^{\beta_q}$ such that $0 \le \beta_q \le \text{val}_q \ell_r$ for all $q \in \text{irred}(\ell_r)$:

(u) let $F$ be the union of the $\mathcal{F}(q)$ over all $q \in \text{irred}(\ell_r)$ such that $\beta_q > 0$,
(v) for $p \in \text{irred}(\ell_0)$, set $a_p$ to 0 if $p \in F$, to $\text{val}_p \ell_0$ otherwise,

(w) for all $A := \prod_{p \in \text{irred}(\ell_0)} p^{\alpha_p}$ such that $0 \le \alpha_q \le a_p$ for all $p \in \text{irred}(\ell_0)$,

(1) compute $\tilde{L}(t, M)$ by (6.8),

...

Here, by the definition of $\mathcal{B}(\pi, q)$, those $q$ for which $\beta_q > 0$ are either in $\pi$ or in $\text{irred}(\ell_r) \setminus \mathcal{R}$. Therefore, the set $F$ computed at step (u) is $\mathcal{F}(\pi)$, because $\mathcal{F}(q) = \emptyset$ if $q \notin \mathcal{R}$. The integer $a_p$ defined at step (v) reflects the definition of $\mathcal{A}(p, \pi)$.

Additionally, the set $\mathcal{F}(q)$ for a given $q$ can be computed efficiently by using the Gr"affe operator $G$ defined by $Gf = \text{Res}_y(y^b - x, f(y))$ (cf. CDDM 2018, §3.2, especially Lemma 3.1), as expressed in the following lemma. The point of this approach is to avoid carrying out divisions of high-degree polynomials $M^s p$ by the polynomials $q$.

**Lemma 7.1.** For any monic irreducible factor $q$ of $\ell_r$, one has

$$ \mathcal{F}(q) = \text{irred}(\ell_0) \cap \{q, \sqrt{G}q, \dots, \sqrt{G}^{r-1}q\}, $$

where $\sqrt{G}f$ denotes the squarefree part of $Gf$, computed by forcing exponents to 1 in the factored representation of $Gf$.

*Proof.* From the definition of $G$ follow the multiplicativity formula $G(fg) = GfGg$ and the degree-preservation formula $\text{deg} Gp = \text{deg} p$. In particular, $G$ preserves divisibility. One has the crucial relations $u | M^i G^i u$ and $G^i M^i u = u^{b^i}$ for any monic irreducible polynomial $u$ and $i \in \mathbb{N}$. In terms of $\sqrt{G}$, the formulas become $u | M^i \sqrt{G}^i u$ and $\sqrt{G}^i M^i u = u$ for $i \in \mathbb{N}$.

Fix an integer $s$ such that $0 \le s < r$. If $q | M^s p$, then upon applying $\sqrt{G}^s$ we find $\sqrt{G}^s q | \sqrt{G}^s M^s p = p$, so $\sqrt{G}^s q = p$. Conversely, if $\sqrt{G}^s q = p$, then upon applying $M^s$ we find $q | M^s \sqrt{G}^s q = M^s p$, so $q | M^s p$. Now, a monic irreducible factor $p$ of $\ell_0$ is forbidden by a monic irreducible factor $q$ of $\ell_r$ if and only if there exists an integer $s$ such that $0 \le s < r$ and $q | M^s p$, that is, if and only if there exists an integer $s$ such that $0 \le s < r$ and $p = \sqrt{G}^s q$. $\square$

## 7.2. Avoiding redundant pairs

Because some of the coprimality conditions defining bounded Gosper-Petkov-šek forms cannot be taken into account in Algorithm 1 until parametrizations of solutions $C$ have been obtained, solutions are naturally found with repetitions. Here, we alleviate this by predicting some of the repetitions from divisibility conditions on $A$ and $B$. By the design of Algorithm 1, a rational solution $u(x)$ will be generated multiple times (before the final step that removes redundancy) if and only if two triples $(A, B, C)$ and $(A', B', C')$ of monic polynomials satisfying

$$ \frac{MC}{C} \frac{M^{r-1}A}{B} = \frac{MC'}{C'} \frac{M^{r-1}A'}{B'} \qquad (7.1) $$
are considered for the same $\zeta$ during the run of the algorithm.

Recall that in Algorithm 1, ramification is taken care of by introducing a new indeterminate $t$ that plays the role of some appropriate root of $x$. The results described in §7.1 for polynomials in $\mathbb{L}[x]$ apply with trivial modifications for polynomials of $\mathbb{L}[t]$, and for the end of §7, we continue with the indeterminate $t$. In particular, we let $\ell_0$ and $\ell_r$ denote $\ell_0(t)$ and $\ell_r(t)$, respectively.

We will study several scenarios in which, if the pair $(A, B)$ is considered in the run of the algorithm and, for a certain $\zeta$, leads to a polynomial $C$ and a rational solution $u$, then another triple $(A', B', C')$ leads to the same $u$ for the same $\zeta$. We will thus obtain rules that attempt to: minimize the multiplicity of $t$ in $A$, maximize the multiplicity of $t$ in $B$, remove factors of the form $Mp/p$ from $A$, introduce factors of the form $Mp/p$ in $B$, replace a divisor $Mp$ with the factor $p$ in $A$, replace the factor $p$ with a divisor $Mp$ in $B$. In all cases, the polynomial $C$ will be adjusted to compensate for the change. This compensation will operate by multiplications/divisions by elements of $\mathbb{L}[t]$ obtained independently of $\zeta$. This makes the actual value of $\zeta$ irrelevant when considering equalities of the form (7.1), which justifies our pruning independently of $\zeta$.

*Pruning rules.* The following four situations provide instances of (7.1) that will permit us to discard pairs $(A, B)$. Note that in each case, monic $A, B, C$, and $p$ induce monic $A', B',$ and $C'$:

* If $B = p\tilde{B}$, define $A' = A$, $B' = Mp\tilde{B}$, $C' = pC$ to get:

$$ \frac{MC}{C} \frac{M^{r-1}A}{B} = \frac{MpMC}{MpC} \frac{M^{r-1}A'}{p\tilde{B}} = \frac{MC'}{C'} \frac{M^{r-1}A'}{B'} $$

* If $A = Mp\tilde{A}$, define $A' = p\tilde{A}$, $B' = B$, $C' = M^{r-1}pC$ to get:

$$ \frac{MC}{C} \frac{M^{r-1}A}{B} = \frac{MC}{C} \frac{M^{r-1}(Mp\tilde{A})}{B} = \frac{M(M^{r-1}pC)M^{r-1}(p\tilde{A})}{M^{r-1}pC} = \frac{MC'}{C'} \frac{M^{r-1}A'}{B'} $$

* If $p \mid Mp$, define $A' = A$, $B' = (Mp/p)B$, $C' = pC$ to get:

$$ \frac{MC}{C} \frac{M^{r-1}A}{B} = \frac{MpMC}{pC} \frac{M^{r-1}A}{(Mp/p)B} = \frac{MC'}{C'} \frac{M^{r-1}A'}{B'} $$

* If $p \mid Mp$ and $A = (Mp/p)\tilde{A}$, define $A' = \tilde{A}$, $B' = B$, $C' = M^{r-1}pC$ to get:

$$ \frac{MC}{C} \frac{M^{r-1}A}{B} = \frac{MC}{C} \frac{M^{r-1}((Mp/p)\tilde{A})}{B} = \frac{M(M^{r-1}pC)M^{r-1}A'}{M^{r-1}pC} = \frac{MC'}{C'} \frac{M^{r-1}A'}{B'} $$
For any monic divisor $A$ of $\ell_0$, any monic divisor $B$ of $\ell_r$, these guarded formulas provide a rule to discard the pair $(A, B)$ if there exists a monic irreducible polynomial $p$ satisfying any of the predicates:

(P1) $B$ is of the form $p\tilde{B}$ and $Mp\tilde{B}$ divides $\ell_r$.

(P2) $A$ is of the form $Mp\tilde{A}$ and $p\tilde{A}$ divides $\ell_0$.

(P3) $p$ divides $Mp$ and $(Mp/p)B$ divides $\ell_r$.

(P4) $p$ divides $Mp$ and $A$ is of the form $(Mp/p)\tilde{A}$.

As an example, we legitimate the use of (P4) after assuming the last situation listed above. Because $A' | A | \ell_0$ and $B' = B | \ell_r$, if $(A, B, C)$ is considered in the algorithm for a certain $\zeta$ and leads to a solution $u$, then the pair $(A', B')$ is also considered, and the solution $u$ found with $(A, B)$ after the polynomial solving step gets $C$ will also be found with $(A', B')$ and the same $\zeta$ after the polynomial solving step gets $C'$. Using the pair $(A', B')$ instead of $(A, B)$ is better, because it potentially leads to more solutions, including polynomials $C'$ not divisible by $M^{r-1}p$. The use of the other predicates is justified by a similar reasoning.

It is clear that the monic nature of $p$ is no restriction: it is needed to relate a monic $B$ and a monic $\tilde{B}$ in (P1), and a monic $A$ and a monic $\tilde{A}$ in (P2); only $Mp/p$ plays a role in (P3) and (P4) and it depends on $\mathbb{L}p$ only.

Furthermore, the reasoning on the predicates does not rely on the irreducibility of $p$: one easily proves that there exists a (monic) irreducible $p$ making a rule apply if and only if there exists a (monic) not necessarily irreducible $p$ making the same (generalized) rule apply. So it is algorithmically permissible to restrict to (monic) irreducible $p$, and doing so avoids having to consider exponentially more factors $p$.

Finally, making $p = t$ in (P3) and (P4) results in slightly more explicit formulations, respectively:

(P5) $\text{val}_t B \le \text{val}_t \ell_r - b + 1$.

(P6) $\text{val}_t A \ge b - 1$.

Note that making $p = t$ in (P1) results in special cases of (P5), and doing so in (P2) results in special cases of (P6).

*Iteration of the rules.* An optimization for factors is obtained from two distinct monic irreducible polynomials $p$ and $q$ satisfying $Mp = pq$: considering the final result of a repeated use of predicates makes it possible to fix certain multiplicities before entering the loops over $A$ and $B$. If for $k \ge 1$, $q^k B$ divides $\ell_r$, an iterated use of (P3) shows that all of $B, \dots, q^{k-1}B$ can be skipped, and we get the predicate:

(P7) $Mp = pq$ for $p \ne q$ and $\text{val}_q B < \text{val}_q \ell_r$.
Similarly, if for $k \ge 1$, $A$ is of the form $q^k \tilde{A}$ and divides $\ell_0$, an iterated use of (P4) shows that all of $q^k \tilde{A}, \dots, q \tilde{A}$ can be skipped, and we get the predicate:

(P8) $Mp = pq$ for $p \neq q$ and $\text{val}_q A > 0$.

Note that $p \neq q$ requires $p \neq t$ and $q \neq t$.

*Implementation of the discarding rule.* The optimization of loop (B) that we introduced in §7.1 is further refined in Algorithm 3: before (u) in the loop over $B$, we insert a step (t) to take (P1) and (P3) into account, and before (1) in the loop over $A$, we insert a step (0) to take (P2) and (P4) into account. Notwithstanding, the specializations (P5) to (P8) are taken care of at step (Z), by restricting the exponents used in $A$ and $B$ to ranges in which the rules (P1) to (P4) cannot apply.

Because (P3) tests those polynomials $Mp/p$ that divide $\ell_r$, a precomputation at step (k) determines those quotients. For any monic irreducible factor $q$ of any such $Mp/p$, $\sqrt{Gq}$ is equal to $p$, thus restricting the search for $p$ to irred($G\ell_r$). A similar discussion applies to (P4) and step (l).

### 7.3. Avoiding redundant computations of Newton polygons

Given an operator $P$, recall that the $\xi_j$ denote the characteristic polynomials of its upper Newton polygon (Definition 6.2), and that $Z(P)$ is the set of all roots $\zeta$ of these characteristic polynomials, according to (6.9). Step (B)(1) of Algorithm 1 computes the characteristic polynomials $\xi_j(X)$ associated with $\tilde{L}$ for every pair $(A, B)$ of divisors of $\ell_0$ and $\ell_r$. It turns out that the roots of the $\xi_j(X)$ do not depend on $(A, B)$. In this section, the degree of a finite Puiseux series is defined as the maximal (rational) exponent appearing with a nonzero coefficient.

**Lemma 7.2.** The set $Z(\tilde{L})$ is equal to $Z(L)$ for every pair $(A, B)$ of polynomials in $\mathbb{L}[t]$. Moreover, the set $\Delta(\tilde{L})$ of degrees of finite Puiseux series solutions of $\tilde{L}y = 0$ is obtained from the set $\Delta(L)$ of degrees of finite Puiseux series solution of $Ly = 0$ by an affine transformation:

$$ \Delta(\tilde{L}) = \mathrm{sh}(\Delta(L)) \quad \text{where} \quad \mathrm{sh}(\delta) = -\frac{b^{r-1} \deg A - \deg B}{b-1} + b^{r-1} \delta. \quad (7.2) $$

*Proof.* The upper Newton polygon of $L$, respectively $\tilde{L}$, is the upper part of the convex hull of the points $(b^k, d_k)$, $0 \le k \le r$, respectively of the points $(b^k, D_k)$, $0 \le k \le r$, with the relation

$$ D_k = b^{r-1}d_k + b^{r-1}\frac{b^k-1}{b-1}\deg A + \frac{b^r-b^k}{b-1}\deg B = b^{r-1}d_k + \alpha b^k + \beta $$

that results directly from (6.8) after setting

$$ \alpha = \frac{b^{r-1} \deg A - \deg B}{b-1} \quad \text{and} \quad \beta = \frac{b^r \deg B - b^{r-1} \deg A}{b-1}. $$
As the coefficient $b^{r-1}$ is positive, the upper convex hull of the points $(b^k, d_k)$, $0 \le k \le r$, is mapped onto the upper convex hull of the points $(b^k, D_k)$, $0 \le k \le r$. As a consequence, the degrees of the finite Puiseux series solutions of $L$ and $\tilde{L}$ are related by (7.2), since they are the opposite of the slopes of the edges (CDDM 2018, Lemma 2.5). Moreover, as the polynomials $A$ and $B$ are monic, the coefficients of the monomials $x^{d_k}M^k$ in $L$ and $x^{D_k}M^k$ in $\tilde{L}$ are the same for $0 \le k \le r$. So the characteristic polynomials are the same for the edges of both upper Newton polygons. $\square$

A first consequence of Lemma 7.2 is that we can change the structure of Algorithm 1: instead of recomputing $Z(\tilde{L})$ at step (2) for each new $(A, B)$, a single computation of $Z(L)$ is done before entering loop (B), and step (3) becomes a loop over $\zeta \in Z(L)$. In the latter loop, we still need to compute a bound on the degree of a solution $C$ for each pair $(A, B)$, but it is deduced from the lemma by considering the largest nonnegative integer in $\Delta(\tilde{L})$ as obtained from $\Delta(L)$ by the parametrization (7.2).

# Part III: Algorithm by reconstruction from series

## 8. Solutions as syzygies

We develop another approach to finding all ramified rational function solutions of the Riccati Mahler equation (R). In comparison to the adaption of Petkovšek's algorithm in §6, which focuses on fixing the possible singularities of a solution of (R), this second approach will proceed in a guess-and-check manner, refining calculations on approximate solutions of the linear equation (L) until a whole set of candidate solutions is proven to be exact solutions. This approach will be embodied in our Algorithm 4.

Our algorithm will search for the set $\mathfrak{R}_{\bar{K}(x^{1/*})}$ of solutions with coefficients in $\bar{K}$, but for the first part of the theory, we more generally consider an intermediate field $\bar{K} \supseteq L \supseteq \bar{K}$, as in §6. We compute Hermite-Padé approximants for auxiliary problems, before recombining them to obtain structured approximants. These auxiliary problems show no logarithm and ramification, making it possible to focus on formal power series calculations. Each auxiliary problem solves a modified Riccati equation for its rational solutions (in $L(x)$) whose series expansions have nonnegative valuation and their leading coefficient is equal to 1.

Let us focus on the space of formal power series solutions of the linear equation (L) and consider an $L$-basis $(z_1, \dots, z_N)$, from which we deduce solutions $u$ of (R) in $L[[x]]$ by the parametrization

$$ u = \frac{M(a_1z_1 + \dots + a_Nz_N)}{a_1z_1 + \dots + a_Nz_N}, \quad a = (a_1 : \dots : a_N) \in \mathbb{P}(\mathbb{L}^N). \quad (8.1) $$

We need to determine those $a$ for which the series $u$ lies in $L(x)$. Suppose that $u$ is indeed some $P/Q$ in $L(x)$, with coprime $P$ and $Q$. After canceling
denominators and using linearity, we obtain the equation

$$(-a_1 P) z_1 + \cdots + (-a_N P) z_N + (a_1 Q) M z_1 + \cdots + (a_N Q) M z_N = 0. \quad (8.2)$$

This linear relation with polynomial coefficients from $\mathbb{L}[x]$ between the series $z_1, \dots, z_N, Mz_1, \dots, Mz_N$ from $\mathbb{L}[[x]]$ is a special case of a linear relation

$$P_1 z_1 + \cdots + P_N z_N + Q_1 M z_1 + \cdots + Q_N M z_N = 0, \quad (8.3)$$

for polynomial coefficients $P_i$ and $Q_i$ from $\mathbb{L}[x]$. A second level of relaxation is obtained, for any $\sigma \in \mathbb{N}$, by Hermite-Padé approximants $(P_1, \dots, P_N, Q_1, \dots, Q_N)$ in $\mathbb{L}[x]^{2N}$, to order $\sigma$, that is, approximate linear relations of the form

$$P_1 z_1 + \cdots + P_N z_N + Q_1 M z_1 + \cdots + Q_N M z_N = O(x^{\sigma}). \quad (8.4)$$

Recombining exact relations (8.3) into structured relations of the form (8.2) reduces to solving a polynomial system in $a_1, \dots, a_N$ (see Definition 8.26), so that we can in principle go back from (8.3) to (8.2). Unfortunately, given bounds on the degrees of the $P_i$ and $Q_i$, it is not clear how to find an accuracy $\sigma$ such that (8.4) implies (8.3). Instead, Algorithm 4 computes the approximate relations (8.4) compatible with our degree bounds on $P$ and $Q$ and attempts to reconstruct hypergeometric solutions starting from these relations. This proceeds by guessing candidate relations for increasing $\sigma$ and rejecting wrong ones until we get no false solutions. We will show (Theorem 9.2) that this process eventually yields all hypergeometric solutions, and nothing but hypergeometric solutions.

The remainder of §8 is dedicated to the details of this strategy, leading to an algorithm presented in §9.1. The discussion involves a number of intermediate objects whose definitions are summarized in Table 1. In terms of the notation in this table, our approach proceeds by computing an overapproximation $\mathcal{V}^{[\sigma]}$ of the limit $\mathcal{A}^{[\infty]}$ we have to compute. At the limit, the relaxed cone $\mathcal{V}^{[\infty]}$ and the characteristic cone $\mathcal{A}^{[\infty]}$ are equal, and are equal to a direct sum of vector spaces. Because, by means of Gröbner basis calculations, we are able to detect nonlinear components of the relaxed cone $\mathcal{V}^{[\sigma]}$, we can develop an algorithm by rejection, which increases $\sigma$ and refines $\mathcal{V}^{[\sigma]}$ until we get the exact set of solutions to (R). A plot summary of the analysis to prove our approach is the following:

* we show that the situation is trivial unless $1 \le \rho^{[\sigma]} \le 2N - 1$;

* we describe the relaxed cone $\mathcal{V}^{[\sigma]}$ as the vanishing set of a polynomial system $\Sigma^{[\sigma]}$ that reflects the degeneracy of the augmented generating matrix $W_a^{[\sigma]}$;

* some case distinction is then needed:

- if $\rho^{[\sigma]} = 2N - 1$, all the series solutions of (L) are hypergeometric and $\mathcal{V}^{[\sigma]} = \mathbb{L}^N$;
<table><caption>Table 1: Inventory of the objects to be introduced in §8 and §9. Each of the quantities X<sup>[σ]</sup> above has a limit when σ goes to ∞, denoted X<sup>[∞]</sup>.</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>σ</td><td></td><td>the truncation order of solution series z<sub>i</sub> to (L) that serve to obtain approximate syzygies to the same order; increased when the algorithm detects it cannot conclude</td></tr><tr><td>Γ<sup>[σ]</sup></td><td>Def. 8.3</td><td>L[x]-module of approximate syzygies (“approximate syzygy module”); nonincreasing and ultimately decreasing w.r.t. σ; constant rank</td></tr><tr><td>Γ<sup>[σ]</sup></td><td>Def. 8.6</td><td>vectorial truncation of the approximate syzygy module to total degree at most B<sub>∞</sub>; nonincreasing and ultimately constant w.r.t. σ; a finite-dimensional L-vector space</td></tr><tr><td>ρ<sup>[σ]</sup></td><td>Def. 8.10</td><td>rank of the module L[x]Γ<sup>[σ]</sup>; nonincreasing and ultimately constant w.r.t. σ</td></tr><tr><td>Δ<sup>[σ]</sup></td><td>Def. 8.13, Prop. 8.19</td><td>L-cone of the a for which structured approximate syzygies exist (“characteristic cone of structured approximate syzygies”, “characteristic cone” for short); nonincreasing and ultimately constant w.r.t. σ; converges to a direct sum of L-vector spaces</td></tr><tr><td>Ψ<sup>[σ]</sup></td><td>Def. 8.20</td><td>L-cone that is a relaxation of Δ<sup>[σ]</sup>, thus generally providing only candidates a (“relaxed cone of structured approximate syzygies”, “relaxed cone” for short); converges to the same limit as Δ<sup>[σ]</sup></td></tr><tr><td>W<sup>[σ]</sup></td><td>Def. 8.6</td><td>matrix computed as a minimal basis from the z<sub>i</sub> (“generating matrix”); its rows generate L[x]Γ<sup>[σ]</sup></td></tr><tr><td>W<sub>a</sub><sup>[σ]</sup></td><td>Def. 8.18</td><td>augmentation of the matrix W<sup>[σ]</sup> (“augmented generating matrix”); singular if and only if a ∈ Υ<sup>[σ]</sup></td></tr><tr><td>Σ<sup>[σ]</sup></td><td>Def. 8.26</td><td>polynomial system obtained from W<sub>a</sub><sup>[σ]</sup>; generates an ideal ⟨Σ<sup>[σ]</sup>⟩ whose variety is Υ<sup>[σ]</sup></td></tr><tr><td>√⟨Σ<sup>[σ]</sup>⟩</td><td>Thm. 9.1</td><td>radical of the ideal ⟨Σ<sup>[σ]</sup>⟩; used to test linearity of the irreducible components of Υ<sup>[σ]</sup></td></tr></tbody></table>

- if ρ<sup>[σ]</sup> < 2N - 1, we isolate the irreducible components of Υ<sup>[σ]</sup>; if one irreducible component is described by a nonlinear system, σ is too small and the algorithm needs to restart with an increased value;

- in both cases, the relaxed cone Υ<sup>[σ]</sup> is described as a union of linear spaces, and there corresponds to each of those spaces a parametrization of a set of candidate solutions to (R);

- each linear parametrization is tested by substituting into the left-hand side of (R): a nonzero evaluation invalidates σ, and the algorithm restarts with an increased value;
• otherwise, the parametrizations thus obtained describe exactly all solutions, up to possible inclusions of a parametrization into another.

*Remark 8.1.* The same idea can be adapted to find the rational solutions of a linear Mahler equation. We further comment on this in §9.2. The latter, as well as §8.7 on the rank of syzygy modules, are only incidental to the flow of the text.

*Remark 8.2.* The problem of determining those $a$ in (8.1) making $u$ rational simplifies a lot when $N=1$, as it reduces to the question whether $Mz_1/z_1$ is rational. Using Hermite-Padé approximation together with the bounds derived in §5 gives an efficient and simple procedure in that context.

## 8.1. Approximate syzygies

In order to describe the linear relations (8.3), hereafter called syzygies, in connection with the approximate linear relations (8.4), hereafter called approximate syzygies, we recall and specialize to our setting some useful notions and facts from the classical theory of Hermite-Padé approximants.

**Definition 8.3.** Given an integer $m \in \mathbb{N}$ and a tuple $f = (f_1, \dots, f_m) \in \mathbb{L}[[x]]^m$, an element $w \in \mathbb{L}[x]^m$ is a syzygy of $f$ when $w \cdot f^T = 0$. Given an integer $\sigma \in \mathbb{N}$, an approximate syzygy of $f$ (to order $\sigma$) is an element $w \in \mathbb{L}[x]^m$ satisfying $w \cdot f^T = O(x^\sigma)$. The module of syzygies is denoted $\mathcal{S}^{[\sigma]}$. The module of approximate syzygies to order $\sigma$ is denoted $\mathcal{S}^{[\sigma]}$.

It is clear that the sequence of the approximate syzygy modules $\mathcal{S}^{[\sigma]}$ is nonincreasing, meaning that for any $\tau \le \sigma$, $\mathcal{S}^{[\tau]} \supseteq \mathcal{S}^{[\sigma]}$. And because a series is 0 if and only if it is $O(x^\sigma)$ for all $\sigma \in \mathbb{N}$, the limit of the sequence is the syzygy module $\mathcal{S}^{[\infty]}$.

*Remark 8.4.* For an approximate syzygy $w$ of $f$ there exists a series $q \in \mathbb{L}[[x]]$ satisfying $(w, q) \cdot (f, x^\sigma)^T = 0$. Therefore, the approximate syzygy module is in general larger than the projection to its first $m$ components of the syzygy module of $(f_1, \dots, f_m, x^\sigma)$, as the latter would restrict the implied $q$ to polynomials.

It is well known (Derksen 1994) that the approximate syzygy module $\mathcal{S}^{[\sigma]}$ is a free $\mathbb{L}[x]$-module of rank $m$, and that this module admits a basis of a specific shape, called a *minimal basis* in the recent literature.

**Definition 8.5.** A basis $(w_1, \dots, w_m)$ of some free $\mathbb{L}[x]$-module is called a minimal basis if for each $i$, the $m$-tuple $w_i = (p_1, \dots, p_m)$ satisfies $\deg p_j \le \deg p_i > \deg p_{j'}$ whenever $j \le i < j'$.

Minimal bases are also well known to be Gröbner bases of the approximate syzygy module for a term-over-position (TOP) ordering (Neiger 2016). Fast algorithms for their calculation have been provided, most notably in (Beckermann and Labahn 1994); see also (Beckermann and Labahn 2000).

To link the previous definitions with the notation of the introduction of §8 and in particular (8.2), (8.3), and (8.4), we set $m := 2N$ and, for $1 \le i \le N$, $f_i := z_i$ and $f_{i+N} := Mz_i$. If the Riccati equation has nontrivial solutions, then
for a large enough $\sigma$ we expect the existence of an approximate syzygy of the form

$$w_0 := (-a_1 P, \dots, -a_N P, a_1 Q, \dots, a_N Q) \qquad (8.5)$$

that is in fact an exact syzygy satisfying $\deg P \le B_{\text{num}}$ and $\deg Q \le B_{\text{den}}$ for the bounds introduced in Proposition 5.2.

Let us consider a minimal basis $(w_1^{[\sigma]}, \dots, w_{2N}^{[\sigma]})$ of the approximate syzygy module $\mathcal{S}^{[\sigma]}$. Any $w \in \mathcal{S}^{[\sigma]}$ reduces to zero by reduction by the minimal basis viewed as a Gröbner basis for the TOP order. Thus, there are polynomials $R_i$ such that

$$w = R_1 w_1^{[\sigma]} + \dots + R_{2N} w_{2N}^{[\sigma]}$$

with $\deg(R_i w_i^{[\sigma]}) \le \deg w$ for $1 \le i \le 2N$: for the quotient $R_i$ to be nonzero, we need $\deg w_i^{[\sigma]} \le \deg w$. In particular, for any $d \in \mathbb{N}$, any $w$ of degree at most $d$ reduces to zero by the elements of the minimal basis themselves of degree at most $d$. With our goal to be able to generate structured syzygies $w_0$ of degree bounded as in (8.5), what just precedes justifies retaining only those elements of the minimal basis whose degrees are at most

$$B_{\infty} := \max(B_{\text{num}}, B_{\text{den}}). \qquad (8.6)$$

**Definition 8.6.** The vectorial truncation of the approximate syzygy module to degree $B_{\infty}$ is the $\mathbb{L}$-vector space

$$\mathcal{T}^{[\sigma]} := \mathcal{S}^{[\sigma]} \cap \mathbb{L}[x]_{\le B_{\infty}}^{2N} = \sum_{\deg w_i^{[\sigma]} \le B_{\infty}} \mathbb{L}[x]_{\le B_{\infty} - \deg w_i^{[\sigma]} w_i^{[\sigma]}}. \qquad (8.7)$$

We write $W^{[\sigma]}$ for the matrix whose rows are those $w_i^{[\sigma]}$ with degree at most $B_{\infty}$, for $1 \le i \le 2N$, and we call it the generating matrix.

Note that the minimal basis remains implicit in the notation of $\mathcal{T}^{[\sigma]}$ and $W^{[\sigma]}$. This is no problem in practice as, once $\sigma$ is fixed, our analysis and algorithm consider a single minimal basis to order $\sigma$.

**Remark 8.7.** Instead of $\mathcal{S}^{[\sigma]} \cap \mathbb{L}[x]_{\le B_{\infty}}^{2N}$ in (8.7), one could suggest using $\mathcal{S}^{[\sigma]} \cap (\mathbb{L}[x]_{\le B_{\text{num}}}^N \times \mathbb{L}[x]_{\le B_{\text{den}}}^N)$, which fits the degree structure of expected exact syzygies (8.5) more tightly. While that should allow for some savings in the initial calculations, the algorithm would anyway have to deal with much higher degrees when the kernel computation (step (C)(h)(6)(ε)[iv] in Algorithm 5) would result in an output of degree bounded by $N \times (B_{\text{num}} + B_{\text{den}}) \le 2N \times B_{\infty}$.

We have already said that, for general $f \in \mathbb{L}[x]^{2N}$, the sequence of the approximate syzygy modules $\mathcal{S}^{[\sigma]}$ is nonincreasing and the rank of $\mathcal{S}^{[\sigma]}$ remains constant at $2N$. We now prove that, for our specific $f$ defined from $z_i$, the sequence is also ultimately decreasing. To this end, without loss of generality we can assume that $v := \text{val} z_1 \le \text{val} z_i$ for $1 \le i \le N$. Thus, for $\sigma \ge v$, $(x^{\sigma-v}, 0, \dots, 0)$ is in $\mathcal{S}^{[\sigma]}$ but not in $\mathcal{S}^{[\sigma+1]}$.
**Lemma 8.8.** The limit $\mathcal{T}^{[\infty]} = \bigcap_{\sigma \in \mathbb{N}} \mathcal{T}^{[\sigma]}$ is equal to $\mathcal{S}^{[\infty]} \cap \mathbb{L}[x]^{2N}_{\le B_{\infty}}$, and the module it spans, $\mathbb{L}[x]\mathcal{T}^{[\infty]}$, is an $\mathbb{L}[x]$-submodule of $\mathcal{S}^{[\infty]}$.

*Proof.* The sequence of the $\mathcal{T}^{[\sigma]}$ is nonincreasing, and, as a subspace of a fixed finite-dimensional vector space, is ultimately constant. The result follows after intersecting with $\mathbb{L}[x]^{2N}_{\le B_{\infty}}$. $\square$

**Remark 8.9.** Because the degree bound $B_{\infty}$ holds a priori only for relations on $(z_i, Mz_i)$, the module $\mathbb{L}[x]\mathcal{T}^{[\infty]}$ can be a strict submodule of the syzygy module $\mathcal{S}^{[\infty]}$.

Altogether, we have just proven the next lemma, which we state after a useful definition.

**Definition 8.10.** For any $\sigma$, possibly $\infty$, let $\rho^{[\sigma]}$ denote the rank of the module $\mathbb{L}[x]\mathcal{T}^{[\sigma]}$.

The rank $\rho^{[\sigma]}$ is also the height of the generating matrix $W^{[\sigma]}$ whose rows in $\mathbb{L}[x]^{2N}$ generate the module (see Definition 8.6).

**Lemma 8.11.** For each $\sigma \in \mathbb{N}$, the $\mathbb{L}[x]$-module $\mathbb{L}[x]\mathcal{T}^{[\sigma]}$ is generated by the rows of the generating matrix $W^{[\sigma]}$. The sequence of these modules is nonincreasing and ultimately constant with limit the submodule $\mathbb{L}[x]\mathcal{T}^{[\infty]}$ of $\mathcal{S}^{[\infty]}$. A similar property holds for the sequence of their ranks, which we write

$$2N \ge \rho^{[\sigma]} \ge \rho^{[\infty]}. \qquad (8.8)$$

**Remark 8.12.** Given a tuple $f = (f_1, \dots, f_m) \in \mathbb{L}[[x]]^m$, some precision $\sigma$, and a minimal basis $w = (w_1, \dots, w_m)$ of the corresponding $\mathbb{L}[x]$-module $\mathcal{S}^{[\sigma]}$ of approximate syzygies, it is immediate, by reduction, that the basis $w$ is also a minimal basis of its $\mathbb{L}'[x]$-module of approximate syzygies at precision $\sigma$ for any superfield $\mathbb{L}'$ of $\mathbb{L}$.

## 8.2. Structured syzygies

Given an operator $L \in \mathbb{K}[x]\langle M \rangle$ of order $r$ with a space of formal power series solutions assumed of $L$-dimension $N > 0$, we fix a basis $z := (z_1, \dots, z_N)$ of the space of series solutions. For any $\sigma \in \mathbb{N}$, the $\mathbb{L}[x]$-module $\mathbb{L}[x]\mathcal{T}^{[\sigma]}$, which consists of unstructured (approximate) syzygies, may contain structured syzygies. We consider in the next definition the set of those $a \in (\mathbb{L}^N)_{\ne 0}$ that correspond to a structured syzygy, which together with 0 make an $L$-cone (see Definition 2.3).

**Definition 8.13.** For each $\sigma \in \mathbb{N}$, the characteristic cone of structured approximate syzygies, or characteristic cone for short, is the $L$-cone

$$\mathcal{A}^{[\sigma]} = \{0\} \cup \{a \in (\mathbb{L}^N)_{\ne 0} \mid \exists P \in \mathbb{L}[x]_{\ne 0}, \exists Q \in \mathbb{L}[x]_{\ne 0}, \\ (-Pa, Qa) \in \mathbb{L}[x]\mathcal{T}^{[\sigma]}\}. \quad (8.9)$$

As a consequence of Lemma 8.11, the sequence of the characteristic cones $\mathcal{A}^{[\sigma]}$ is nonincreasing and ultimately constant, with limit a cone $\mathcal{A}^{[\infty]}$ that we proceed to describe.
**Lemma 8.14.** Let $s_1, s_2, \dots$ denote the dimensions of the nontrivial similarity classes $(\mathfrak{H}_1)_{\neq 0}, (\mathfrak{H}_2)_{\neq 0}, \dots$ of the hypergeometric series solutions of $\mathcal{L}$. There exist full-rank matrices $H_1, H_2, \dots$ such that:

* the $\mathcal{L}$-vector spaces $\mathcal{L}^{s_i} H_i$ are in direct sum in $\mathcal{L}^N$;

* the limiting characteristic cone $\mathcal{A}^{[\infty]}$ is the union of the vector spaces $\mathcal{L}^{s_i} H_i$.

*Proof.* Consider an $\mathcal{L}(x)$-similarity class of hypergeometric series solutions $\mathfrak{H}_{\neq 0}$ of $(\mathcal{L})$, described by a basis $h := (h_1, \dots, h_s)$. Expressing this family in the basis $z$ provides us with a matrix $H \in \mathcal{L}^{s \times N}$ such that $h^T = \mathcal{H}z^T$. Each $c \in (\mathcal{L}^s)_{\neq 0}$ yields a rational-function solution $M(ch^T)/(ch^T)$ of $(\mathbb{R})$, also equal to $M(c\mathcal{H}z^T)/(c\mathcal{H}z^T)$. After writing $P/Q$ for this rational function, we get that $(-PcH, QcH)$ is a structured approximate syzygy to any order $\sigma$, so that $cH$ is in the cone $\mathcal{A}^{[\infty]}$. This cone therefore contains $\mathcal{L}^s H$. Iterating over the nontrivial similarity classes $(\mathfrak{H}_1)_{\neq 0}, (\mathfrak{H}_2)_{\neq 0}, \dots$ of hypergeometric series solutions, we therefore get dimensions $s_1, s_2, \dots$ and matrices $H_1, H_2, \dots$ such that each $\mathcal{L}^{s_i} H_i$ is contained in $\mathcal{A}^{[\infty]}$. By point 1 of Theorem 2.9 applied to $D = \mathfrak{D}_{\mathcal{L}}$ (see Definition 3.1) and $F = \mathcal{L}((x))$, the $\mathfrak{H}_i$ are in direct sum in $\mathfrak{D}_{\mathcal{L}}$, inducing, by the map $a \mapsto az^T$, that the $\mathcal{L}^{s_i} H_i$ are in direct sum in the $\mathcal{L}$-space generated by $\mathcal{A}^{[\infty]}$. Conversely, any nonzero $a \in \mathcal{A}^{[\infty]}$ yields a nonzero rational solution $P/Q$ of the Riccati Mahler equation, thus, by point 2 of Theorem 2.9 must come from some $(\mathfrak{H}_i)_{\neq 0}$. So, the cone $\mathcal{A}^{[\infty]}$ consists solely of the union of the $\mathcal{L}$-vector spaces $\mathcal{L}^{s_i} H_i$. $\square$

In view of (8.8), if some rank $\rho^{[\sigma]}$ is 0, then $\rho^{[\infty]} = 0$, meaning there are no nontrivial syzygies at all, and our algorithm will terminate in this case. So, we continue the analysis by assuming $\rho^{[\infty]} > 0$. Additionally, if the rank $\rho^{[\infty]}$ were $2N$, there would exist some finite $\sigma$ for which $\rho^{[\sigma]} = 2N$, which would yield a nonsingular square matrix $W^{[\sigma]}$. We could then perform row operations over $\mathcal{L}[x]$ so as to produce $2N$ nonzero polynomials $p_i$ satisfying $p_i z_i = p_{N+i} M z_i = 0$ for $1 \le i \le N$, but this would be a contradiction to the non-nullity of the $z_i$. Thus, we have $0 < \rho^{[\infty]} < 2N$. Consequently, we also have $\rho^{[\sigma]} < 2N$ for large enough $\sigma$. If, conversely, we observe $\rho^{[\sigma]} = 2N$, then we know that $\sigma$ is too small to provide a description of the wanted solutions, and our algorithm will continue with increased $\sigma$.

**Example 8.15.** We use as a working example throughout §8 the Mahler operator of order $r = 4$ and degree $d = 258$ introduced in Example 1.5. This operator admits a four-dimensional vector space of formal power series solutions in $\mathbb{Q}[[x]]$. For any $\sigma \le 120$, the computation of the generating matrix $W^{[\sigma]}$ (Definition 8.6) leads to a rank $\rho^{[\sigma]} = 8 = 2N$ (Definition 8.10), proving that such $\sigma$ are too low to provide candidates.

**Hypothesis 8.16.** For §8.3, §8.4, and §8.5, we assume $1 \le \rho^{[\infty]}$.

### 8.3. A linear-algebra condition on structured syzygies

From now on, we assume $1 \le \rho^{[\sigma]} \le 2N - 1$. The following lemma is a variant of Cramer's rule.
**Lemma 8.17.** Let $M$ be a matrix of size $(n+1) \times n$ and rank $n$, with coefficients in a field $F$. Then, the left kernel of $M$ has dimension 1 over $F$ and a nonzero kernel element is $K := (\Delta_1, \dots, (-1)^{i+1}\Delta_i, \dots, (-1)^n\Delta_{n+1})$, where $\Delta_i$ denotes the determinant of the square submatrix obtained by removing the $i$th row.

*Proof.* Augmenting $M$ by any column $C$ of it on its left yields a matrix with determinant 0. Expanding the determinant with respect to the first column shows that $K$ satisfies $KC = 0$. Combining all columns $C$ yields $KM = 0$. Because rk $M = n$, the minors $\Delta_i$ cannot all be 0 simultaneously, which implies $K \neq 0$. $\square$

The characteristic cone $\mathcal{A}^{[\sigma]}$ was defined by (8.9). We proceed to formulate an alternative description of it in terms of the generators of the module $\mathbb{L}[x]\mathcal{T}^{[\sigma]}$ provided by the rows of $W^{[\sigma]}$.

**Definition 8.18.** Let $W_a^{[\sigma]}$ denote the matrix obtained by augmenting $W^{[\sigma]}$ at its bottom with the two-row matrix

$$ \begin{pmatrix} a & 0 \\ 0 & a \end{pmatrix} = \begin{pmatrix} a_1 & \dots & a_N & 0 & \dots & 0 \\ 0 & \dots & 0 & a_1 & \dots & a_N \end{pmatrix}. \qquad (8.10) $$

We call it the augmented generating matrix.

**Proposition 8.19.** The characteristic cone $\mathcal{A}^{[\sigma]}$ (Definition 8.13) satisfies

$$ \mathcal{A}^{[\sigma]} = \{0\} \cup \{a \in (\mathbb{L}^N)_{\neq 0} \mid \exists P \in \mathbb{L}[x]_{\neq 0}, \exists Q \in \mathbb{L}[x]_{\neq 0}, \exists \Pi \in \mathbb{L}[x]^{\rho^{[\sigma]}}, (\Pi_1, \dots, \Pi_{\rho^{[\sigma]}}, P, -Q) W_a^{[\sigma]} = 0\}. \quad (8.11) $$

Note that this characterization is in fact independent of the choice of the generating matrix $W^{[\sigma]}$. We introduce the cone $\mathcal{Y}^{[\sigma]}$ defined analogously by the looser constraint $(\Pi, P, -Q) \neq 0$ in place of $P \neq 0$ and $Q \neq 0$: this is the cone of values of $a$ for which $W_a^{[\sigma]}$ has a nontrivial left kernel.

**Definition 8.20.** For each $\sigma \in \mathbb{N}$, the relaxed cone of structured approximate syzygies, or relaxed cone for short, is the $\mathbb{L}$-cone

$$ \mathcal{Y}^{[\sigma]} = \{0\} \cup \{a \in (\mathbb{L}^N)_{\neq 0} \mid \exists K \in (\mathbb{L}[x]^{\rho^{[\sigma]}+2})_{\neq 0}, KW_a^{[\sigma]} = 0\}. \quad (8.12) $$

Note the inclusion $\mathcal{A}^{[\sigma]} \subseteq \mathcal{Y}^{[\sigma]}$ of the characteristic cone into the relaxed cone. The following lemma describes a case of equality.

**Lemma 8.21.** For $\sigma$ large enough, meaning $\mathcal{T}^{[\sigma]} = \mathcal{T}^{[\infty]}$, assume $1 \le \rho^{[\sigma]} \le 2N - 1$. Then:

1. if, for some nonzero $a$, the row $K = (\Pi_1, \dots, \Pi_{\rho^{[\sigma]}}, -P, Q) \neq 0$ satisfies $KW_a^{[\sigma]} = 0$, then neither $P$ nor $Q$ can be zero;

2. $\mathcal{A}^{[\sigma]} = \mathcal{Y}^{[\sigma]} = \mathcal{A}^{[\infty]}$.
**Proof.** Taking $a$ as in the first point implies $Paz^T = QaMz^T + \Pi W^{[\sigma]}(z, Mz)^T$, which reduces to $Paz^T = QaMz^T$ because $W^{[\sigma]}(z, Mz)^T = 0$ by the hypothesis on $\sigma$. Now, if $P$ is zero, then $Qa = 0$ because the entries of $Mz$ are independent over $\mathbb{L}$, forcing $Q = 0$ because $a \neq 0$. Similarly, $Q = 0$ implies $P = 0$. So, if $P$ or $Q$ is zero, then both must be zero, which contradicts that the generating matrix $W^{[\sigma]}$ has independent rows over $\mathbb{L}(x)$, thus proving $\mathcal{V}^{[\sigma]} \subseteq \mathcal{A}^{[\sigma]}$, and therefore $\mathcal{A}^{[\sigma]} = \mathcal{V}^{[\sigma]}$. Moreover, if two integers $\sigma_1$ and $\sigma_2$ are such that $\mathcal{T}^{[\sigma_1]} = \mathcal{T}^{[\sigma_2]}$, then $\mathcal{A}^{[\sigma_1]} = \mathcal{A}^{[\sigma_2]}$ by the definition (8.9). In particular, for all $\sigma$ satisfying $\mathcal{T}^{[\sigma]} = \mathcal{T}^{[\infty]}$, all the cones $\mathcal{A}^{[\sigma]}$ are equal and equal to $\mathcal{A}^{[\infty]}$. $\square$

**Lemma 8.22.** The limit $\mathcal{V}^{[\infty]}$ of relaxed cones exists and is equal to the limit $\mathcal{A}^{[\infty]}$ of characteristic cones.

**Proof.** The existence of $K$ in the definition of $\mathcal{V}^{[\sigma]}$ is equivalent to the existence of a nonzero element in

$$ \mathbb{L}[x]^{\rho^{[\sigma]}} W^{[\sigma]} \cap \mathbb{L}[x]^2 \begin{pmatrix} a & 0 \\ 0 & a \end{pmatrix} = \mathbb{L}[x] \mathcal{T}^{[\sigma]} \cap \mathbb{L}[x]^2 \begin{pmatrix} a & 0 \\ 0 & a \end{pmatrix}. $$

By Lemma 8.11, this intersection is nonincreasing and ultimately constant, so that $\mathcal{V}^{[\sigma]}$ is nonincreasing and ultimately constant as well. This shows the existence of $\mathcal{V}^{[\infty]}$. Now, by Hypothesis 8.16, if $\rho^{[\sigma]} < 2N$ for some $\sigma$, then for all $\tau \ge \sigma$, we have $1 \le \rho^{[\tau]} \le 2N - 1$. So if some $\sigma$ satisfies the hypotheses of Lemma 8.21, then all $\tau \ge \sigma$ satisfy them as well. We deduce $\mathcal{V}^{[\infty]} = \mathcal{A}^{[\infty]}$. $\square$

**Remark 8.23.** In order to compute $\mathcal{A}^{[\infty]}$, Algorithm 4 will not manipulate $\mathcal{A}^{[\sigma]}$ directly and will manipulate $\mathcal{V}^{[\sigma]}$ instead. It will therefore work by increasing $\sigma$ until the equalities of Lemma 8.21(2) are satisfied.

## 8.4. Candidate solutions

Two cases need different analyses: (i) $\rho^{[\sigma]} = 2N - 1$; (ii) $\rho^{[\sigma]} \le 2N - 2$.

### 8.4.1. Case $\rho^{[\sigma]} = 2N - 1$.

In the first case, the $(\rho^{[\sigma]} + 2) \times (\rho^{[\sigma]} + 1)$-matrix $W_a^{[\sigma]}$ has a nontrivial left kernel for any value of $a$: the relaxed cone $\mathcal{V}^{[\sigma]}$ is the whole space $\mathbb{L}^N$. Additionally, the $\rho^{[\sigma]} \times (\rho^{[\sigma]} + 1)$-matrix $W^{[\sigma]}$ has full rank over $\mathbb{L}(x)$, and thus has a 1-dimensional right kernel. Applying Lemma 8.17 to the transposed matrix provides a polynomial generator $K^T = (\Delta_1, \dots, (-1)^{i+1}\Delta_i, \dots, -\Delta_{2N})^T$ of this right kernel, given by minors $\Delta_i$ obtained by removing columns. As ranks are unchanged by extending the base field to $\mathbb{L}((x))$, the generating matrix $W^{[\sigma]}$ also has a 1-dimensional right kernel over $\mathbb{L}((x))$, and the same $K^T$ is a generator of this $\mathbb{L}((x))$-vector space.

**Lemma 8.24.** For $\sigma$ large enough, meaning $\mathcal{T}^{[\sigma]} = \mathcal{T}^{[\infty]}$, assume $\rho^{[\sigma]} = 2N - 1$. Then:

$$ \mathcal{A}^{[\sigma]} = \mathcal{V}^{[\sigma]} = \mathbb{L}^N, $$
* those rational solutions $u \in \mathbb{L}(x)$ to the Riccati equation (R) that are quotients of the form $Mw/w$ for $w = \sum_{i=1}^{N} a_i z_i$ and $a \in (\mathbb{L}^N)_{\neq 0}$ are given by the parametrization

$$ (a_1 : \dots : a_N) \mapsto u = \frac{\sum_{i=1}^{N} a_i (-1)^{i+N+1} \Delta_{i+N}}{\sum_{i=1}^{N} a_i (-1)^{i+1} \Delta_i}, \quad \text{with } a \in \mathbb{P}^{N-1}(\mathbb{L}). \qquad (8.13) $$

**Proof.** By the assumption on $\sigma$, $(z_1, \dots, z_N, Mz_1, \dots, Mz_N)^T$ must be in the kernel of the generating matrix $W^{[\sigma]}$, and thus must be a multiple of $K^T$ by a series $z_0 \in \mathbb{L}((x))_{\neq 0}$. In particular, $z_i = (-1)^{i+1} \Delta_i z_0$ for $1 \le i \le N$, so all the $z_i$ are $\mathbb{L}(x)$-similar. For any $a \in (\mathbb{L}^N)_{\neq 0}$, the nonzero series $w = \sum_{i=1}^{N} a_i z_i$ is a solution to (L) and an immediate calculation shows that the quotient $Mw/w$ is a rational solution from $\mathbb{L}(x)$ given by the formula in (8.13). That is, $w$ is hypergeometric and, by point 2 of Theorem 2.9, $u := Mw/w$ is a solution to (R). This shows that (8.13) parametrizes all solutions $u$ that are quotients $Mw/w$ and rational. Any such solution, given by $a \neq 0$, provides nonzero $P$ and $Q$ satisfying $Mw/w = P/Q$. The nullity of $(P, -Q)(\begin{smallmatrix} a & 0 \\ 0 & a \end{smallmatrix})(z, Mz)^T$ implies that $(P, -Q)(\begin{smallmatrix} a & 0 \\ 0 & a \end{smallmatrix})$ is a linear combination of the rows of $W^{[\sigma]}$ with rational function coefficients, that is, there exists $\Pi \in \mathbb{L}(x)^{\rho^{[\sigma]}}$ satisfying $(\Pi, P, -Q)W_a^{[\sigma]} = 0$. After clearing denominators, this shows $a \in \mathcal{A}^{[\sigma]}$. So we proved the equality $\mathcal{A}^{[\sigma]} = \mathbb{L}^N$, from which derives the equality $\mathcal{Y}^{[\sigma]} = \mathbb{L}^N$ because $\mathcal{A}^{[\sigma]} \subseteq \mathcal{Y}^{[\sigma]} \subseteq \mathbb{L}^N$. $\square$

The algorithm does not decide at first if $\mathcal{T}^{[\sigma]} = \mathcal{T}^{[\infty]}$. Instead, the formula (8.13) of Lemma 8.24 introduces a parametrized candidate for rational solutions. If this parametrized rational function is later verified to be an actual parametrized solution, in other words if the algorithm reaches a $\sigma$ large enough to ensure

$$ u = \frac{\sum_{i} a_{i} z_{i+N}}{\sum_{i} a_{i} z_{i}}, $$

then the similarity class of hypergeometric solutions of the $z_i$ is proved to have $\mathbb{L}$-dimension $N$.

**Example 8.25** (continuing from Example 8.15). For each $\sigma$ between 121 and 123, the computation of the generating matrix $W^{[\sigma]}$ leads to a rank $\rho^{[\sigma]} = 7 = 2N - 1$. In each case, the rational candidate (8.13) has a numerator of degree 81 in $x$ and a denominator of degree 83, both involving linear parameters $a_1, \dots, a_4$. It will turn out, however, that this candidate need not parametrize only solutions (see Example 8.30).

8.4.2. Case $1 \le \rho^{[\sigma]} \le 2N - 2$.

In the second case, the matrix $W_a^{[\sigma]}$ is of size $(\rho^{[\sigma]} + 2) \times (2N)$ with $\rho^{[\sigma]} + 2 \le 2N$. This matrix is full rank over $\mathbb{L}(x, a_1, \dots, a_N)$ and we want to find all choices of $a = (a_1, \dots, a_N)$ in $\mathbb{L}^N$ that make the rank drop. They are described by canceling all maximal minors. These minors are polynomials in $a$ and $x$ that are homogeneous of degree 2 in the $a_i$. We gather their coefficients with respect to $x$ into a system $\Sigma^{[\sigma]}$ of polynomials of $\mathbb{L}[a]$.
**Definition 8.26.** Let $\Sigma^{[\sigma]} \subseteq \mathbb{L}[a]$ denote the polynomial system consisting of the coefficients according to the monomial basis $(x^m)_{m \in \mathbb{N}}$ of the maximal minors of the augmented generating matrix $W_a$.

So, the cone $\mathcal{V}^{[\sigma]}$ of all $a$ that make $W_a^{[\sigma]}$ rank deficient is the variety of common zeros in $\mathbb{L}^N$ of $\Sigma^{[\sigma]}$. The following lemma makes the description of the equality case $\mathcal{A}^{[\sigma]} = \mathcal{V}^{[\sigma]}$ in Lemma 8.21 more explicit.

**Lemma 8.27.** For $\sigma$ large enough, meaning $\mathcal{T}^{[\sigma]} = \mathcal{T}^{[\infty]}$, assume $1 \le \rho^{[\sigma]} \le 2N - 2$. Then, for any nonzero $a$ of $\mathbb{L}^N$ that makes $W_a^{[\sigma]}$ rank deficient, the rank over $\mathbb{L}(x)$ of $W_a^{[\sigma]}$ is $\rho^{[\sigma]} + 1$.

*Proof.* Because rk $W^{[\sigma]} = \rho^{[\sigma]}$, the rank rk $W_a^{[\sigma]}$ is either $\rho^{[\sigma]}$ or $\rho^{[\sigma]} + 1$. In the former case, the first $\rho^{[\sigma]}$ rows of $W_a^{[\sigma]}$ generate its image $\mathbb{L}(x)^{\rho^{[\sigma]}+2}W_a^{[\sigma]}$, and there exists $K = (\Pi_1, \dots, \Pi_{\rho^{[\sigma]}}, -P, 0) \in \mathbb{L}[x]^{\rho^{[\sigma]}+2}$ with $P \ne 0$ satisfying $KW_a^{[\sigma]} = 0$. Consequently, $Paz^T = (\Pi_1 w_1^{[\sigma]} + \dots + \Pi_{\rho^{[\sigma]}} w_{\rho^{[\sigma]}})(z, Mz)^T$, which is 0 because $W^{[\sigma]}(z, Mz)^T = 0$ by the hypothesis on $\sigma$. It follows that $az^T = 0$, which contradicts that $z$ is an $\mathbb{L}$-basis. So the rank must be $\rho^{[\sigma]} + 1$. $\square$

By Lemma 8.14, the limit cone $\mathcal{A}^{[\infty]}$ is a union of linear spaces. The relaxed cone $\mathcal{V}^{[\sigma]}$ approximates the characteristic cone $\mathcal{A}^{[\sigma]}$ by containing it and is equal to it for sufficiently large $\sigma$. Our algorithm therefore computes a primary decomposition of the radical of the ideal generated by $\Sigma^{[\sigma]}$, then tests whether each prime ideal thus obtained describes a linear space.

There exist algorithms for computing primary decompositions over an algebraically closed field like $\mathbb{K}$. Algorithms to be found in the literature (Gianni, Trager, and Zacharias 1988; Decker, Greuel, and Pfister 1999) return primary ideals represented by (minimal reduced) Gröbner bases, and we will prove by Theorem 9.1 that testing linearity of the irreducible components of the relaxed cone $\mathcal{V}^{[\sigma]}$ amounts to checking that all returned ideals are generated by linear forms.

So, on the one hand, if any element of any of the (minimal reduced) Gröbner bases is nonlinear, then our algorithm restarts with a larger $\sigma$.

On the other hand, if all Gröbner bases define linear spaces, our algorithm extracts a more explicit parametrization of candidates $u = P/Q$. To this end, take each Gröbner basis in turn, solve the system it represents for $a$ to determine parameters $(g_1, \dots, g_v)$ and some matrix $S$ of size $v \times N$ such that $a = gS$. The following lemma shows the resulting form of candidate solutions $P/Q$, parametrized by $g$.

**Lemma 8.28.** Assume $1 \le \rho^{[\sigma]} \le 2N - 2$ and let $g \in \mathbb{L}^v$ satisfy $gS \in (\mathcal{V}^{[\sigma]})_{\ne 0}$. Then, the left kernel of $W_a^{[\sigma]}$ is generated over $\mathbb{L}(x)$ by $K = (\Pi_1, \dots, \Pi_{\rho^{[\sigma]}}, -P, Q)$ where $P$ and $Q$ are nonzero and depend linearly on $g$, and the $\Pi_i$ are quadratic in $g$.

*Proof.* The matrix $W_a^{[\sigma]}$ has size $(\rho^{[\sigma]} + 2) \times (2N)$ and by Lemma 8.27 it has $\mathbb{L}(x)$-rank $\rho^{[\sigma]} + 1$. As a consequence, we can select $\rho^{[\sigma]} + 1$ linearly independent
columns to obtain a matrix $\Omega(g)$ of rank $\rho^{[\sigma]} + 1$, made of an upper $\rho^{[\sigma]} \times (\rho^{[\sigma]} + 1)$ block extracted from $W_a^{[\sigma]}$ and independent of $g$, and of two rows depending linearly on $g$. Lemma 8.17 provides a nonzero $K = (\Pi_1, \dots, \Pi_{\rho^{[\sigma]}}, -P, Q)$ satisfying $K\Omega(g) = 0$, with nonzero $P$ and $Q$ by Lemma 8.21. Since the columns of $\Omega(g)$ generate those of $W_a^{[\sigma]}$, the product $KW_a^{[\sigma]}$ is zero as well. The definition of $K$ by minors and the degree structure of $W_a^{[\sigma]}$ with respect to $a$ provide the degrees in $g$ of its entries. $\square$

**Example 8.29** (continuing from Example 8.25). For $\sigma = 124$, computing the generating matrix $W^{[\sigma]}$ leads to a rank $\rho^{[\sigma]} = 6 = 2N - 2$. Augmenting $W^{[\sigma]}$ by stacking the two-row matrix (8.10) and taking the coefficients with respect to $x$ of the single minor to be considered yields a system $\Sigma^{[\sigma]}$ of 30 polynomials of degree 2 in $a_1, \dots, a_4$. The corresponding irredundant prime decomposition is $\mathcal{J}_1 \cap \mathcal{J}_2 \cap \mathcal{J}_3$ for $\mathcal{J}_1 = \langle a_1 - a_3, a_2 - a_4 \rangle$, $\mathcal{J}_2 = \langle a_1 - a_2, a_1 + a_3, a_1 + a_4 \rangle$, $\mathcal{J}_3 = \langle a_1 - a_4, a_1 + a_2, a_1 + a_3 \rangle$. The associated linear subspaces of $\mathbb{L}^4$, predicted to be in direct sum by Lemma 8.14, have respective dimensions 2, 1, and 1. In the present case, these dimensions add up to $N = 4$. The ideals $\mathcal{J}_2$ and $\mathcal{J}_3$ yield isolated candidates $P/Q$ (projective dimension 0), while the first ideal yields a candidate parametrized by a projective line (projective dimension 1). Focusing on $\mathcal{J}_1$, we solve it as $(a_1, a_2, a_3, a_4) = (g_1, g_2) \begin{pmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \end{pmatrix}$. The specialized augmented matrix becomes

$$ W_a^{[\sigma]} = \begin{pmatrix} x & -1 & x & -1 & 0 & 0 & 0 & 0 \\ p_{2,1}^{(36)} & p_{2,2}^{(37)} & p_{2,3}^{(0)} & p_{2,4}^{(37)} & p_{2,5}^{(0)} & p_{2,6}^{(1)} & p_{2,7}^{(0)} & p_{2,8}^{(1)} \\ 0 & 0 & -1 & 0 & x^2 & 0 & 1 & x \\ 0 & 0 & 0 & -1 & 0 & x^2 & x & 1 \\ -1 & 0 & 0 & 0 & 1 & x & x^2 & 0 \\ 0 & -1 & 0 & 0 & x & 1 & 0 & x^2 \\ g_1 & g_2 & g_1 & g_2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & g_1 & g_2 & g_1 & g_2 \end{pmatrix}, $$

where the polynomials $p_{2,j}^{(d_j)}$ in the second row have degrees $d_j$ for $1 \le j \le 8$. This matrix has $\rho^{[\sigma]} + 2 = 8$ rows and we observe that the first 7 columns are linearly independent over $\mathbb{L}(x, g_1, g_2)$, leading to an $8 \times 7$ matrix $\Omega(g)$ that plays the same role as the matrix with same name in the proof of Lemma 8.28. The determinants of the submatrices of $\Omega(g)$ obtained by deleting the $j$th row for $1 \le j \le 8$ are all of the form $\Delta_j = q\delta_j$ with $q$ a degree 36 polynomial and

$$ \begin{align*} \delta_1 &= -(g_1^2 - g_2^2)x, & \delta_2 &= 0, \\ \delta_3 &= \delta_5 = -(g_1 - g_2x + g_1x^2)(g_1 + g_2x), & \delta_4 &= \delta_6 = (g_2 - g_1x + g_2x^2)(g_1 + g_2x), \\ \delta_7 &= -(g_1 + g_2x^3), & \delta_8 &= -(1 + x^2 + x^4)(g_1 + g_2x). \end{align*} $$

The left kernel is the $\mathbb{L}(x)$-line generated by $K := (\Delta_1, \dots, \Delta_7, -\Delta_8)$. This
provides a rational candidate

$$ \frac{P}{Q} = \frac{\Delta_7}{\Delta_8} = \frac{g_1 + g_2 x^3}{g_1 + g_2 x} \frac{1}{1 + x^2 + x^4}. \qquad (8.14) $$

### 8.5. Validating candidate solutions

Whether obtained when $\rho^{[\sigma]} = 2N - 1$ (§8.4.1) or $\rho^{[\sigma]} \le 2N - 2$ (§8.4.2), a candidate $u = P/Q$ has its numerator and denominator parametrized linearly by some $g = (g_1, \dots, g_v)$ (where $g = a$ and $v = N$ if in the former case). So a candidate is really a parametrization by $g$ in some $\mathbb{L}^{v \ne 0}$ of a family of rational functions from $\mathbb{L}(x)$. Such a candidate parametrization is obtained algorithmically as a rational function over a transcendental extension of $\mathbb{L}$ generated by symbolic generators that represent the parameters $g$.

Suppose a candidate parametrization could provide genuine solutions for some specializations of $g$ as well as rational functions that do not solve (R) for some other specializations of $g$. For sure, enlarging $\sigma$ would ultimately reject the wrong specializations. However, one could try to avoid these new calculations and hope to separate the correct parameters from the wrong ones by using just the data at hand. For example, one could try to substitute $P/Q$ for $u$ in (R) and identify coefficients to 0, but this would lead to a polynomial system of degree $r$ in $g$, which we would not know how to solve efficiently.

So, verifying that a sufficiently large $\sigma$ has been used amounts to verify for each parametrized candidate that $P/Q$ is truly a solution of (R) for all choices of $g$. As a first criterion, any fixed parametrized candidate in normal form $P/Q$ can thus be rejected if either $\deg P > B_{\text{num}}$ or $\deg Q > B_{\text{den}}$, for the bounds defined by (5.1) and (5.2). Otherwise, after substituting $P/Q$ for $u$ in the left-hand side of (R), either we obtain 0 and $P/Q$ is a valid solution, or $\sigma$ has to be increased.

**Example 8.30** (continuing from Example 8.29). Gathering the results from the previous examples, we see that the rank $\rho^{[\sigma]}$ of the submodule $\mathbb{L}[x] \mathcal{F}^{[\sigma]}$ defined by (8.7) of the approximate syzygy module $\mathcal{S}^{[\sigma]}$ is nonincreasing when $\sigma$ increases. Meanwhile, the cone $\mathcal{A}^{[\sigma]}$ of parameters of candidate solutions decreases as well, until it reaches the limit $\mathcal{A}^{[\infty]}$ that corresponds to true solutions only. Let us describe this in more detail:

* $\rho^{[\sigma]} = 2N = 8$ for $\sigma \le 120$. We already explained that these $\sigma$ are too small.

* $\rho^{[\sigma]} = 2N - 1 = 7$ for $121 \le \sigma \le 123$. As the degrees in $x$ of each candidate are $81 > 344/9 = B_{\text{num}} \simeq 38.2$ for the numerator and $83 > B_{\text{den}} = 86/3 \simeq 28.7$ for the denominator, we reject those candidates.

* $\rho^{[\sigma]} = 6 \le 2N - 2$ for $124 \le \sigma \le 126$. For $\sigma = 124$ and the parametrized candidate (8.14), substituting into the Riccati equation proves that (8.14) defines a family of true rational solutions parametrized by $(g_1 : g_2) \in \mathbb{P}^1(\mathbb{Q})$. Similarly, the ideals $\mathfrak{I}_2$ and $\mathfrak{I}_3$ of Example 8.29 each yields one solution. As all candidates are solutions, $\mathcal{A}^{[\sigma]} = \mathcal{A}^{[\infty]}$ for this $\sigma$, and therefore for all $\sigma \ge 124$.
* $\rho^{[\sigma]} = 5 \le 2N - 2$ for $127 \le \sigma$. As we will see in §8.7, each similarity class of hypergeometric solutions of dimension $s$ (over $\mathbb{L} = \mathbb{Q}$) contributes $2s - 1$ linearly independent syzygies (of degree at most $B_{\infty} = B_{\text{num}} \simeq 38.2$), and in presence of several classes, the $2s - 1$ add up to a lower bound on the rank $\rho^{[\infty]}$ of the submodule $\mathbb{L}[x]\mathcal{T}^{[\infty]}$ of the module of syzygies on hypergeometric solutions. In the present example, this lower bound on $\rho^{[\infty]}$ is $(2 \times 2 - 1) + (2 \times 1 - 1) + (2 \times 1 - 1) = 5$. A computation for $\sigma = 127$ finds $\rho^{[\sigma]} = 5$, proving that $\rho^{[\sigma]} = \rho^{[\infty]} = 5$, as well as $\mathcal{T}^{[\sigma]} = \mathcal{T}^{[\infty]}$, for all $\sigma \ge 127$.

The reader might be surprised to note that the equality $\mathcal{A}^{[\sigma]} = \mathcal{A}^{[\infty]}$ is satisfied before $\mathcal{T}^{[\sigma]} = \mathcal{T}^{[\infty]}$ occurs: indeed, the definition of $\rho^{[\sigma]}$ as the rank of $\mathbb{L}[x]\mathcal{T}^{[\sigma]}$ and the inequality $\rho^{[124]} = 6 > \rho^{[127]} = 5 = \rho^{[\infty]}$ imply that $\mathcal{T}^{[124]}$ is strictly larger than $\mathcal{T}^{[\infty]}$, whereas $\mathcal{A}^{[124]} = \mathcal{A}^{[\infty]}$.

## 8.6. Ramified rational solutions

As a consequence of Theorem 3.22 (and Definition 3.21), all the ramified rational solutions of (R) with a given leading coefficient $\lambda$ are images by $y \mapsto My/y$ of the $\mathbb{L}$-vector space of solutions $y$ of (L) in $(\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/q_\lambda}))$ for the ramification bound $q_\lambda$ given by (3.10), and this set is nonempty if and only if $\lambda \in \mathbb{L} \cap \Lambda'$. When nonempty, this is $\mathfrak{R}_{\mathbb{L},\lambda}$.

Let us fix a $\lambda \in \mathbb{L} \cap \Lambda'$. Then, we first change the operator $L = L(x, M)$ into $L(x, \lambda M)$ to reduce the search in $(\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/q_\lambda}))$ to a search in $\mathbb{L}((x^{1/q_\lambda}))$, and we next apply the transformation explained in (CDDM 2018, §2.7, particularly Lemma 2.21) to reduce the search in $\mathbb{L}((x^{1/q_\lambda}))$ to a search in $\mathbb{L}[x]$. To make the operator resulting from those combined transformations more explicit: consider the rightmost $\lambda$-admissible edge of the lower Newton polygon of $L$ whose denominator is coprime with $b$ (see Definition 3.6); express its slope as $-p_\lambda/q_\lambda$, which is possible by the definition of $q_\lambda$; and denote by $c_\lambda$ its intercept on the ordinate axis. By defining

$$L_{\lambda}(x, M) = x^{-q_{\lambda}c_{\lambda}} L(x^{q_{\lambda}}, \lambda M) x^{p_{\lambda}}, \quad (8.15)$$

where the factor $x^{-q_\lambda c_\lambda}$ simply ensures that the coefficients of $L_\lambda$ are co-prime polynomials, we obtain a linear Mahler operator whose solutions $z \in \mathbb{L}[x]$ parametrize the solutions $y$ of $L$ in $(\ln x)^{\log_b \lambda} \mathbb{L}((x^{1/q_\lambda}))$ by $y(x) = (\ln x)^{\log_b \lambda} x^{p_\lambda/q_\lambda} z(x^{1/q_\lambda})$. This parametrization is bijective, owing to (CDDM 2018, Prop. 2.19). To obtain all the ramified rational solutions $u$ of (R), it is therefore sufficient to apply the method of the previous subsections to each $L_\lambda$ and to obtain $u = My/y$ for $y$ defined from $z$ as above.

**Remark 8.31.** An important fact, to be used when applying Theorem 9.1 in the proof of Theorem 9.2, is that the solution space of $L_\lambda$ in $\mathbb{L}((x^{1/*}))$ is in fact included in $\mathbb{L}[x]$.

**Remark 8.32.** The set $\Lambda$ is clearly computable (see (3.7)). The previous reasoning shows that its subset $\Lambda'$ (Definition 3.21) is also computable: determining the series solutions of $L_\lambda$ for all $\lambda \in \Lambda$ will decide which of the $\lambda$ is in $\Lambda'$. However,
the calculation of $\Lambda'$ will remain implicit in Algorithm 4: it will loop over $\lambda \in \Lambda$ and just continue with the next $\lambda$ if no series solution of $L_\lambda$ is found.

**Example 8.33.** To show an example of calculation of $L_\lambda$, let us consider the radix $b \ge 2$ and the linear Mahler operator $L = \lambda_0 x^\omega - M + M^r$, with $\lambda_0$ an algebraic number, $\omega$ a positive integer coprime to $b-1$, and $r$ an integer larger than 1. The characteristic polynomial of the leftmost edge of the lower Newton polygon (Definition 3.6) has a single root, $\lambda_0$, and this edge is the only $\lambda_0$-admissible edge. The change described by (8.15) results in the linear operator $L_{\lambda_0} = \lambda_0(1 - M + \lambda_0^{r-1}x^{(b^r-b)\omega}M^r)$. The power series solutions of $L_{\lambda_0}$ are multiples of

$$z = 1 - \lambda_0^{r-1} x^{(b^r-b)\omega} - \lambda_0^{r-1} x^{(b^r-b)b\omega} + \dots, \quad (8.16)$$

from which we deduce a line of solutions of $L$ in $\mathcal{D}$ generated by

$$y = (\ln x)^{\log_b \lambda_0} x^{\omega/(b-1)} \left( 1 - \lambda_0^{r-1} x^{(b^r-b)\omega/(b-1)} + \dots \right).$$

### 8.7. A supplementary remark on the rank of syzygy module

Strikingly, the module of syzygies contains but is not limited to the direct sum over similarity classes of the module of syzygies contributed by each class. This section describes the potential interactions between the latter. It is not used in the rest of the article.

The following lemma quantifies the syzygies that appear in each similarity class.

**Lemma 8.34.** For an $\mathbb{L}(x)$-similarity class $\mathfrak{H} \neq \mathfrak{H}_0$ of $\mathbb{L}(x)$-hypergeometric elements, let $(y_1, \dots, y_s)$ denote a basis of the $\mathbb{L}$-vector space $\mathfrak{H}$ and $S$ the column vector $(y_1, \dots, y_s, My_1, \dots, My_s)^T$. The module of the row vectors $R \in \mathbb{L}[x]^{2s}$ such that $RS = 0$ has rank $2s-1$.

*Proof.* Making $F = \mathbb{L}(x)$ in Lemma 2.8 implies that $\mathfrak{H}$ is an $\mathbb{L}$-vector space. Since the $y_i$, $1 \le i \le s$, lie in the same $\mathbb{L}(x)$-similarity class, we get equalities $y_i = q_i y_1$ for some $q_i$ in $\mathbb{L}(x)$ and all $i$. As $y_1$ is $\mathbb{L}(x)$-hypergeometric, there exists $u_1$ in $\mathbb{L}(x)$ such that $My_1 = u_1 y_1$, from which it results $My_i = (Mq_i)u_1 y_1$. Therefore, $y_1, \dots, y_s, My_1, \dots, My_s$ are $2s$ elements of the line $\mathbb{L}(x)y_1$. The module of their syzygies thus has rank $2s-1$. $\square$

However, other syzygies may be produced by the interaction between several similarity classes. We cannot quantify the phenomenon, and merely give an example.

**Example 8.35.** Let us consider, with $b = 2$, the lshm (least common left multiple) in $\mathbb{Q}[x] \langle M \rangle$ of the operators $L_1 = (1 - 2x^2)M - (1 - 2x)$ and $L_2 = (1 - 3x^2)M - (1 - 3x)$, which respectively annihilate the rational functions $y_1 = 1/(1 - 2x)$ and $y_2 = 1/(1 - 3x)$. This is the operator

$$P_1 := (6x^8 - 5x^4 + 1) M^2 - (6x^6 + 6x^5 + x^4 + -x^3 - 4x^2 + x + 1) M \\ + (6x^4 + x^3 - 4x^2 + x).$$
We then slightly modify it by truncating its coefficients, to obtain the operator

$$P_2 := (-5x^4 + 1) M^2 - (x^4 - 5x^3 - 4x^2 + x + 1) M + (x^3 - 4x^2 + x),$$

in such a way that $P_2$ has a 2-dimensional space of power series solutions, like the operator $P_1$ we started with. One of them is the infinite product

$$y_3 = \prod_{k \ge 0} M^k \frac{1 - 5x^2}{1 - 4x + x^2}$$

and another linearly independent one is

$$y_4 = x + 5x^2 + 19x^3 + 71x^4 + 265x^5 + 983x^6 + 3667x^7 + 13661x^8 + O(x^9).$$

At this point, we can observe that the lclm $L$ of $P_1$ and $P_2$ admits three right factors of order 1, namely $M - u_i$, $1 \le i \le 3$, with

$$u_1 = \frac{1 - 2x}{1 - 2x^2}, \quad u_2 = \frac{1 - 3x}{1 - 3x^2}, \quad u_3 = \frac{1 - 4x + x^2}{1 - 5x^2}.$$

As $y_1$ and $y_2$ are in the same $\mathbb{L}(x)$-similarity class, but not $y_3$, we could (wrongly) expect a dimension $(2 \cdot 2 - 1) + (2 \cdot 1 - 1) = 4$ for the rank of the syzygy module of the column vector $(y_1, \dots, y_4, My_1, \dots, My_4)^T$. However the operator $P_2$ factors as

$$P_2 = (1 - 5x^4)(M - v)(M - u_3), \quad \text{with } v = x \frac{1 - 5x^2}{1 - 5x^4}$$

and the operator $M - v$ admits the rational solution $z = x/(1 - 5x^2)$. As a consequence the series $y_4$, which satisfies $P_2y_4 = 0$, also satisfies the equation $(M - u_3)y_4 = cz$ for some constant $c$. But $z$, as a rational function, lies in the same similarity class as $y_1$ and $y_2$. So we have an additional syzygy, between $y_1$, $y_4$, and $My_4$, and the rank of the syzygy module proves to be 5.

# 9. Algorithms by Hermite-Padé approximation

## 9.1. Algorithm by sieving candidates for solving the Riccati equation

By composing the results developed in §8, we derive Algorithm 4, in which the body of the loop over the approximation order $\sigma$ has been isolated as Algorithm 5. The termination and correctness of the algorithm will be proved in Theorem 9.2, after having proved a structural property of the ideals it manipulates in Theorem 9.1. Before this, we comment on choices and algorithms used at a few specific steps:

* For each $\lambda$, the algorithm computes a basis of power series solutions truncated to an order $\sigma_0$ that is sufficient to distinguish the basis elements. This order is determined at step (C)(c). To explain the calculation, recall
Input: A Riccati Mahler equation (R) with coefficients in $K[x]$.

Output: The set of ramified rational functions $u \in K(x^{1/*})$ that solve (R).

(A) Compute the lower Newton polygon $\mathcal{N}$ of $L := \sum_{k=0}^{r} \ell_k(x)M^k$, then the set $\Lambda$ by (3.8).

(B) Determine the leftmost edge of $\mathcal{N}$ and compute its slope $-\nu$ and intercept $\mu$.

(C) For each $\lambda$ in $\Lambda$,

(a) compute the ramification bound $q_\lambda$ as defined by (3.10),

(b) determine the rightmost $\lambda$-admissible edge of $\mathcal{N}$ and compute its slope $-p_\lambda/q_\lambda$ and intercept $c_\lambda$,

(c) compute the operator $L_\lambda$ by (8.15), the rationals $\nu_\lambda$ and $\mu_\lambda$ given by (9.1) in terms of $\mu$ and $\nu$, and the integer $\sigma_0 := \lfloor \nu_\lambda \rfloor + 1$,

(d) compute a basis $(z_1, \dots, z_N)$ of solutions in $K[\lambda][\lfloor x \rfloor]$ to the equation $L_\lambda z = 0$, truncated to order $O(x^{\sigma_0})$,

(e) if $N=0$, then continue to the next $\lambda$,

(f) compute bounds $B_{\text{num}}$, $B_{\text{den}}$, $B_\infty$ by (5.1), (5.2), and (8.6) applied to $L_\lambda$,

(g) initialize the solution set by $R_\lambda := \emptyset$,

(h) for $\sigma := (\frac{1+\sqrt{5}}{2})^k \sigma_0$ given by successive $k=0, 1, 2, \dots$, execute Algorithm 5 in the current context,

(i) for each $u$ in $R_\lambda$, change $u(x)$ into $x^{(b-1)p_\lambda/q_\lambda}u(x^{1/q_\lambda})$.

(D) Return the union of the sets $R_\lambda$ over $\lambda \in \Lambda$.

**Algorithm 4:** Ramified rational solutions to a Riccati Mahler equation by Hermite-Padé approximants and prime decompositions.

(CDDM 2018, Prop. 2.6 and Lemma 2.21) that the computation of a basis of solutions for the equation $L_\lambda z = 0$ depends on the rational numbers

$$ \nu_{\lambda} = q_{\lambda}\nu - p_{\lambda}, \quad \mu_{\lambda} = q_{\lambda}(\mu - c_{\lambda}), \qquad (9.1) $$

where $-\nu$ and $\mu$ are respectively the slope and intercept of the leftmost edge of the lower Newton polygon of $L$, and, as was explained before (8.15), $-p_\lambda/q_\lambda$ and $c_\lambda$ are respectively the slope and intercept of the rightmost 1-admissible edge of the Newton polygon of $L_\lambda$. The order $\sigma_0$ is then defined to be $\lfloor \nu_\lambda \rfloor + 1$.

* Solving for truncated solutions $z_i$ at step (C)(d) is done by solving a linear system of size $(\lfloor \nu_\lambda \rfloor + 1) \times (\lfloor \mu_\lambda \rfloor + 1)$ as per (CDDM 2018, §2.6, Algorithm 5).

* The initial value of $\sigma$ chosen at step (C)(h) uses all data already available at this point. Next, choosing a geometric sequence to grow $\sigma$, as opposed to, for example, an arithmetic sequence, is justified by the superlinear complexity with respect to $\sigma$ of the calculations in the loop body. The
(1) extend the basis elements $z_1, \dots, z_N$ to solutions truncated to order $O(x^\sigma)$ by
unrolling recurrences,

(2) compute a minimal basis of the module of approximate syzygies of $f = (z_1, \dots, z_N, Mz_1, \dots, Mz_N)^T$ to order $O(x^\sigma)$,

(3) build a matrix $W$ of dimension $\rho \times 2N$ by extracting from the minimal basis
the (independent) rows of degree at most $B_\infty$,

(4) if $\rho = 0$, continue to the next $\lambda$,

(5) if $\rho = 2N - 1$:

(α) for $i = 1, \dots, 2N$, compute the minor $\Delta_i$ obtained after removing the
$ith$ column from $W$,

(β) define $\mathcal{C} := \{u\}$ for the candidate $u$ provided by (8.13), with parameters $a_i$ replaced with $g_i$,

(6) if $1 \le \rho \le 2N - 2$:

(α) let $W_a$ be the $(\rho+2) \times 2N$ matrix obtained by appending the two-row
matrix $\begin{pmatrix} a_1 & \cdots & a_N \\ 0 & \cdots & 0 \\ 0 & \cdots & a_1 \\ \vdots & & \vdots \\ 0 & \cdots & a_N \end{pmatrix}$ below $W$,

(β) let $\langle \Sigma \rangle$ be the ideal of $\overline{\mathbb{K}}[a_1, \dots, a_N]$ generated by the set $\Sigma \subseteq \overline{\mathbb{K}}[\lambda][a_1, \dots, a_N]$ of the coefficients with respect to $x$ of the $\begin{pmatrix} 2N \\ \rho+2 \end{pmatrix}$ minors of order $\rho+2$ of $W_a$,

(γ) compute an irredundant prime decomposition $\bigcap_{j=1}^{s} \mathfrak{p}_j$ of $\sqrt{\langle \Sigma \rangle}$, given
for each $j$ by a Gröbner basis $(p_{j,1}, \dots, p_{j,m(j)})$ of $\mathfrak{p}_j$,

(δ) initialize the candidate set by $\mathcal{C} := \emptyset$,

(ε) for $j$ from 1 to $s$ do

[i] if for some $k$, the polynomial $p_{j,k}$ is nonlinear, continue to the next $\sigma$,

[ii] solve the linear system $\{p_{j,k} = 0\}$ for the unknowns $a_1, \dots, a_N$, so
as to get a full rank matrix $S$ of dimension $v \times N$, with $0 \le v \le N$, and a parametrization $(a_1, \dots, a_N) = (g_1, \dots, g_v)S$,

[iii] if $v=0$, continue to the next $j$,

[iv] compute the left kernel of the matrix $W_q S$,

[v] if the kernel has dimension 2, or if it has dimension 1 and is
generated by a row $K$ with $K_{\rho+1}K_{\rho+2} = 0$, continue to the next $\sigma$,

[vi] add the normalized form of $u := -K_{\rho+1}/K_{\rho+2}$ to $\mathcal{C}$,

(7) for each $u = P/Q$ in $\mathcal{C} \setminus R_\lambda$ do

(α) if $\deg P > B_{\text{num}}$ or if $\deg Q > B_{\text{den}}$, then continue to the next $\sigma$,

(β) if $P/Q$ cancels the left-hand side of the Riccati equation (R), then add $u$ to $R_\lambda$, else continue to the next $\sigma$,

(8) continue to step (C)(i) in Algorithm 4, thus quitting the loop over $\sigma$.

Algorithm 5: Body of step (C)(h) in Algorithm 4.
golden ratio used to increase $\sigma$ is rather arbitrary, but reduces a majority of timings for the set of equations tested.

* Prolonging the truncated solutions $z_i$ at step (C)(h)(1) is done by (CDDM 2018, §2.3, Algorithm 3).

* Computing the minimal basis at step (C)(h)(2) can be done by either the algorithm in (Derksen 1994) or the algorithm in (Beckermann and Labahn 1994). In both cases, the theoretical description admits algebraic coefficients in $\mathbb{K}$.

* At step (C)(h)(6)(γ), the algorithm uses the classical notion of an irreducible primary decomposition of an ideal, which applies in particular to radical ideals, in which case the primary ideals are even prime (Zariski and Samuel 1958, p. 209, especially Theorem 5); we will speak of an irreducible primary decomposition in such a situation. The computation of such a prime decomposition for radical ideals in $\mathbb{K}[a_1, \dots, a_N]$ can be done by a variant of the algorithm in (Gianni, Trager, and Zacharias 1988); see also (Becker and Weispfenning 1993, Algorithms RADICAL, p. 394, and PRIMDEC, p. 396).

The proof of Theorem 9.1 is based on the Nullstellensatz. Consequently, for the rest of §9.1 we specialize $\mathbb{L}$ to $\mathbb{K}$ in the theory developed in §8.

**Theorem 9.1.** Consider an operator $L \in \mathbb{K}[x] \langle M \rangle$ and a basis $z := (z_1, \dots, z_N)$ of its space of formal power series solutions in $\mathbb{K}[[x]]$. Assume that $L$ has no solutions in $\mathbb{K}((x^{1/*}))$ that are not in $\mathbb{K}[[x]]$. Consider as well the family of modules $\mathbb{K}[x] \mathcal{T}^{[\sigma]}$ of ranks $\rho^{[\sigma]}$ and the family of characteristic cones $\mathcal{A}^{[\sigma]}$ defined in Definitions 8.6 and 8.13 (see also Proposition 8.19). Assume that there exists $\sigma$ satisfying the hypothesis $1 \le \rho^{[\sigma]} \le 2N - 2$ of §8.4.2 and large enough to have the relation $\mathcal{T}^{[\sigma]} = \mathcal{T}^{[\infty]}$. The system $\Sigma := \Sigma^{[\sigma]} \subseteq \mathbb{K}[a] := \mathbb{K}[a_1, \dots, a_N]$ of Definition 8.26 generates an ideal $\langle \Sigma \rangle$ whose radical $\sqrt{\langle \Sigma \rangle}$ is a finite intersection of prime ideals, each given by linear polynomials with coefficients in $\mathbb{K}$.

*Proof.* Lemma 8.14 proves the existence of a finite family of ideals $q_i \subseteq \mathbb{K}[a]$ generated by $\mathbb{K}$-linear polynomials, and satisfying

$$ \mathcal{A}^{[\infty]} = \bigcup_{i \in I} V(q_i), \qquad (9.2) $$

where $V(\cdot)$ denotes the variety in $\mathbb{K}^N$ of an ideal. Every $\mathbb{K}$-subspace $\mathfrak{H}$ of $\mathbb{K}((x^{1/*}))$ associated with a $\mathbb{K}(x^{1/*})$-similarity class $\mathfrak{H} \neq 0$ of hypergeometric Puiseux series solutions of $L$ is, by the assumption on the solutions of $L$, in fact included in $\mathbb{K}[[x]]$, and therefore admits a finite basis $h = (h_1, \dots, h_s)$ with elements in $\mathbb{K}((x^{1/*}))$ that are $\mathbb{K}$-linear combinations of the series $z_i$. Therefore, there is a matrix $H \in \mathbb{K}^{s \times N}$ such that $h^T = H z^T$.

Each $V(q_i)$ is in bijection with such an $\mathfrak{H}$ by the $\mathbb{K}$-linear map $a \mapsto az^T$, and more specifically $V(q_i) = \mathbb{K}^s H$. Consequently, $q_i$ is generated by a system of $N - s$ linearly independent linear polynomials with coefficients in $\mathbb{K}$.
We now have the successive equalities

$$V(\sqrt{\langle\Sigma\rangle}) = V(\langle\Sigma\rangle) = \mathcal{Y}^{[\sigma]} = \mathcal{A}^{[\infty]} = \bigcup_{i \in I} V(q_i) = V(\bigcap_{i \in I} q_i),$$

where the first equality if by (Cox, Little, and O’Shea 2015, Theorem 7(ii), p. 183), the second equality is by definition, the third equality results by Lemma 8.21 from the assumption $\mathcal{F}^{[\sigma]} = \mathcal{F}^{[\infty]}$, the fourth equality is (9.2), and the fifth equality is by (Cox, Little, and O’Shea 2015, Theorem 15, p. 196). Retaining the equality between the first and last terms, then passing to ideals, Hilbert’s Nullstellensatz (Cox, Little, and O’Shea 2015, Theorem 2, p. 179) and the definition of a radical provide the first equality in

$$\sqrt{\langle\Sigma\rangle} = \sqrt{\bigcap_{i \in I} q_i} = \bigcap_{i \in I} q_i,$$

where, additionally, the second is because the ideals $q_i$ are prime, therefore radical, and because an intersection of radical ideals is radical ((Zariski and Samuel 1958, Theorem 9, p. 147) or (Cox, Little, and O’Shea 2015, Proposition 16, p. 197)). Finally, any irredundant primary decomposition of $\sqrt{\langle\Sigma\rangle}$ is obtained by retaining a subfamily of the $q_i$, because the latter are prime. This ends the proof as the $q_i$ are defined by linear polynomials. □

**Theorem 9.2.** *Algorithm 4 terminates and correctly computes all solutions of (R) in $\overline{\mathbb{K}}(x^{1/*})$.*

*Proof.* The general structure of the algorithm is a loop over $\lambda$ at step (C), with independent calculations for different $\lambda$, so it is sufficient to prove that for each $\lambda \in \Lambda$ the calculation terminates and computes all solutions $u$ with leading coefficient $\lambda$.

For each $\lambda$, the algorithm introduces the operator $L_\lambda \in \overline{\mathbb{K}}[\lambda][x]\langle M \rangle$ at step (C)(c) and a basis $(z_1, \dots, z_N)$ of its solutions in $\overline{\mathbb{K}}[\lambda][\langle x \rangle]$ at step (C)(d). Note that by Lemma 3.9, the dimension of formal power series solutions of $L_\lambda$ in $\overline{\mathbb{K}}[\langle x \rangle]$ is also equal to $N$, and that the $\overline{\mathbb{K}}[\lambda]$-basis $z$ of solutions in $\overline{\mathbb{K}}[\lambda][\langle x \rangle]$ is also a $\overline{\mathbb{K}}$-basis of the solutions in $\overline{\mathbb{K}}[\langle x \rangle]$. It is therefore sufficient to prove that step (C)(h) computes all rational solutions that can be written $M(az^T)/(az^T)$ for $a \in (\overline{\mathbb{K}}^N)_{\neq 0}$ to get that the algorithm determines all ramified rational solutions $u \in \overline{\mathbb{K}}(x^{1/*})$ of (R) with leading coefficient $\lambda$ at step (C)(i).

With respect to termination, step (C) may quit early at step (C)(e), but if the run goes beyond step (C)(e), then the basis $z$ contains nonzero entries. By the comment for the case $\rho^{[\infty]} = 2N$ after Lemma 8.14, this forces $\rho^{[\infty]} < 2N$ and the algorithm continues with an unbounded inner loop over $\sigma$ at step (C)(h). This inner loop can quit early at step (C)(h)(4), ending the calculation for the current $\lambda$. The inner loop can only be relaunched from step (C)(h)(6)($\varepsilon$)[i] if a nonlinear polynomial if detected, from step (C)(h)(6)($\varepsilon$)[v] if a nonzero rational function cannot be defined, and from steps (C)(h)(7)($\alpha$) and ($\beta$) if a false solution is detected.
At the construction of $z$ at step (C)(d), no $z_i$ is in $O(x^{\sigma_0})$, and therefore no $z_i$ is ever in $O(x^{\sigma})$ after the extension step (C)(h)(1). Consequently, $\rho$ is never $2N$ in the inner loop, which is why only the cases $0 \le \rho \le 2N - 1$ are considered in the algorithm.

The objects computed inside the loop at step (C)(h), most of which have coefficients in $\mathbb{K}[\lambda]$, relate to the theoretical objects defined in §8 for $\mathbb{L} = \bar{\mathbb{K}}$. We already proved that $z$ is a basis for the solutions of $L_\lambda$ in $\bar{\mathbb{K}}[[x]]$. At step (C)(h)(2), the computation of the minimal basis, by Derksen's algorithm (1994), Beckermann and Labahn's algorithm (Beckermann and Labahn 1994) or other known algorithms, is independent of any extension of $\mathbb{K}[\lambda]$ in which the coefficients of the input would be seen. Consequently, the matrix so obtained provides a minimal basis of the approximate syzygy module both over $\mathbb{K}[\lambda][x]$ and over $\bar{\mathbb{K}}[x]$. By the second interpretation, the matrix $W$ computed at step (C)(h)(3) is the matrix $W^{[\sigma]}$ of Definition 8.6 for $\mathbb{L} = \bar{\mathbb{K}}$, and its height, $\rho$ in the algorithm description, is the rank $\rho^{[\sigma]}$ of the theory, also for $\mathbb{L} = \bar{\mathbb{K}}$. For the rest of the proof, we consider the objects $\mathcal{T}^{[\sigma]}$ (vectorial truncation of the approximate syzygy module), $\mathcal{A}^{[\sigma]}$ (characteristic cone), $\mathcal{V}^{[\sigma]}$ (relaxed cone), $\Sigma^{[\sigma]}$, which are all obtained only from $W$, and the theory applied to $\mathbb{L} = \bar{\mathbb{K}}$ (see Definitions 8.6, 8.13, and 8.26, and Proposition 8.19).

Suppose that the algorithm runs without exiting the inner loop, thus making $\sigma$ grow indefinitely. In particular, $\rho$ is never 0, which would cause early quitting, so we indefinitely have $1 \le \rho \le 2N - 1$. From some point on, we have $\mathcal{T}^{[\sigma]} = \mathcal{F}^{[\infty]}$, so that Lemma 8.21 applies, and $\mathcal{V}^{[\sigma]} = \mathcal{A}^{[\sigma]} = \mathcal{A}^{[\infty]}$ holds. If the calculation enters step (C)(h)(5), it continues to step (C)(h)(7). Otherwise, the calculation enters step (C)(h)(6). Then $\rho \le 2N - 2$ and the ideal $\sqrt{\langle\Sigma^{[\sigma]}\rangle}$ computed at step (C)(h)(6)(γ) is the ideal of the variety $\mathcal{V}^{[\sigma]}$. At this point, it is legitimate to apply Theorem 9.1 to the operator $L_\lambda$, which obviously has coefficients in $\bar{\mathbb{K}}$, all other hypotheses being fulfilled in terms of objects over $\mathbb{L} = \bar{\mathbb{K}}$ (see in particular Remark 8.31). We conclude that the prime ideals $\mathfrak{p}_j$ have linear generators, so that the algorithm passes beyond step (C)(h)(6)(ε)[i]. By Lemma 8.27, the kernel computed at step (C)(h)(6)(ε)[iv] has dimension 1, and by Lemma 8.21(1), all $K$ computed satisfy $K_{\rho+1}K_{\rho+2} \ne 0$, so that the algorithm also passes beyond step (C)(h)(6)(ε)[v]. In this case again, the calculation continues to step (C)(h)(7). At this point, because $\mathcal{V}^{[\sigma]} = \mathcal{A}^{[\infty]}$ can only describe true solutions, owing to Lemma 8.14, the validation step (C)(h)(7) finds no false solution and the loop ends, a contradiction.

We have just seen that the loop over $\sigma$ terminates, and that it does so after proving that all candidates are true solutions. Consider the final value of $\sigma$ and exclude the trivial case of early termination at step (C)(h)(4), that is, assume $\rho^{[\sigma]} > 0$. In the inner loop body, the algorithm runs either step (C)(h)(5) or (6). In the first case, $\mathcal{V}^{[\sigma]} = \bar{\mathbb{K}}^N$, which, in view of Lemma 8.24, leads to a rational candidate $u$ parametrized by the free parameter $(g_1, \dots, g_N)$ at step (C)(h)(5)(β). In the second case, the prime decomposition computed at step (C)(h)(6)(γ) represents $\sqrt{\langle\Sigma^{[\sigma]}\rangle}$, so the union of the varieties of the $\mathfrak{p}_j$ in $\bar{\mathbb{K}}^N$ is the relaxed cone $\mathcal{V}^{[\sigma]}$, and the rational candidates are the fractions $P/Q$ implied
by Lemma 8.28. By our hypothesis of a terminating $\sigma$, all the candidate rational solutions obtained previously are verified to be true solutions by step (C)(h)(7), that is $\mathcal{Y}^{[\sigma]} \subseteq \mathcal{A}^{[\infty]}$.

Because of the inclusions $\mathcal{A}^{[\infty]} \subseteq \mathcal{A}^{[\sigma]} \subseteq \mathcal{Y}^{[\sigma]}$ that are valid independently of $\sigma$, we have thus proved $\mathcal{Y}^{[\sigma]} = \mathcal{A}^{[\infty]}$, that is, step (C)(h)(8) returns exactly the set of all solutions $u$ with leading coefficient $\lambda$. $\square$

*Remark 9.3.* Although solving by computing a prime decomposition may exhibit a bad worst-case behavior in theory, we have not investigated how to improve the theoretical complexity of the step after observing experimentally that this is never the bottleneck of execution: the computation of Hermite–Padé approximants always takes more time. In addition, the special structure of a union of spaces makes it plausible that one could develop a better method than the general primary decomposition algorithm.

**Example 9.4.** Let us return to Example 8.33 to show that the number of steps before obtaining the final value of $\sigma$ can be made arbitrarily large in Algorithm 4. The algorithm first determines the parameters $\nu = \omega/(b-1)$ and $\mu = \omega b/(b-1)$. For $\lambda = \lambda_0$, it then gets the ramification bound $q_{\lambda_0} = b-1$, and the rightmost $\lambda_0$-admissible edge (see (3.10) and Definition 3.6) happens to be the leftmost, with slope $-\nu$, thus providing $p_{\lambda_0} = \omega$ and $c_{\lambda_0} = \mu$. This implies $\nu_{\lambda_0} = \mu_{\lambda_0} = 0$. The order $\sigma = \sigma_0$ defined at step (c) for the series expansion at the first iteration of the loop (h) is $\omega + 1$, and solving $L_{\lambda_0}$ yields $N = 1$ at step (d), with a basis $z$ of the form $z = (z_1)$ for $z_1 = 1 + O(x^{\sigma_0})$. The algorithm verifies that the candidate obtained from $z_1 = 1$ is a false solution. In view of (8.16), to distinguish $z_1$ from 1, the algorithm will need to increase the order to approximately $(b^r - b)\omega$. For the sake of exposition, consider a variant setting of the algorithm that multiplies $\sigma$ by $b$ at each iteration (instead of the golden ratio), so that $\sigma = b^k \sigma_0$ (starting with $k=0$). The algorithm then needs to reach $k=r$ before $z$ can be distinguished from 1. For $k=r$, if $\sigma - e > B_{\text{num}}$ for the exponent $e = (b^r - b)\omega$ appearing in (8.16), then the rank $\rho^{[\sigma]}$ is zero and the algorithm immediately stops; if $\sigma - e \ge B_{\text{num}}$, the algorithm continues to increase $k$ until it decides to stop. In all cases, the final value of $\sigma$ is at least $b^r(\omega + 1)$. For a numerical example, set $b=4$, $L = x^{10} - M + M^4$. It takes five iterations and computation with accuracy $O(x^{2816})$ to conclude that the Riccati equation has no solution: this is a case where the algorithm stops at $k=r$, with $2816 = 4^4(10+1) > B_{\text{num}} = 157$ (for the bound $B_{\text{num}}$ relative to $L_{\lambda_0}$).

## 9.2. Rational solving of the linear Mahler equation

In (CDDM 2018), we developed an algorithm to compute the rational solutions of the linear Mahler equation (L). There, we derived degree bounds for the numerators and denominators of rational solutions (see Propositions 3.16, 3.17, and 3.21 in that reference, or, for alternative bounds, (Bell and Coons
2017)). Here<fnref n="4" /> we sketch how the approach by Hermite-Padé approximants of the present §8 adapts, based on these bounds, to give an alternative algorithm for computing rational solutions in the linear case. This adaptation is not necessary for the rest of the article.

To study the linear equation, we modify the definition of the characteristic cone $\mathcal{A}^{[\sigma]}$ and that of the augmented matrix $W_a^{[\sigma]}$ (see Definition 8.18 and Proposition 8.19). We start with the vector of series $(z_1, \dots, z_N, 1)^T = (z, 1)^T$ and, for any given $\sigma$, we compute a minimal basis from which we extract the rows compatible with the degree bounds. This yields a $\rho^{[\sigma]} \times (N+1)$-matrix $W^{[\sigma]}$ satisfying

$$W^{[\sigma]}(z, 1)^T = O(x^{\sigma}).$$

We define $W_a^{[\sigma]}$ as the matrix $W^{[\sigma]}$ augmented at its bottom with the two-row matrix

$$\begin{pmatrix} a & 0 \\ 0 & 1 \end{pmatrix} = \begin{pmatrix} a_1 & \dots & a_N & 0 \\ 0 & \dots & 0 & 1 \end{pmatrix}, \qquad (9.3)$$

and we let

$$\mathcal{A}^{[\sigma]} = \{0\} \cup \{a \in (\mathbb{L}^N)_{\neq 0} \mid \exists P \in \mathbb{L}[x]_{\neq 0}, \exists Q \in \mathbb{L}[x]_{\neq 0}, \exists \Pi \in \mathbb{L}[x]^{\rho^{[\sigma]}}, (\Pi_1, \dots, \Pi_{\rho^{[\sigma]}}, P, -Q)W_a^{[\sigma]} = 0\}. \quad (9.4)$$

The motivation is that for $a \in \mathcal{A}^{[\sigma]}$ and corresponding $(\Pi, P, Q)$, if $W^{[\sigma]}(z, 1)^T$ is zero, we get $(\Pi_1, \dots, \Pi_{\rho^{[\sigma]}}, P, -Q)W_a^{[\sigma]}(z, 1)^T = 0$, implying $Paz^T - Q = 0$, that is, $Q/P = az^T$.

In this context, the rank of the generating matrix $W^{[\sigma]}$ satisfies $\rho^{[\sigma]} \le N$. Again, we distinguish two cases:

* If $\rho^{[\sigma]} = N$, we expect, provided $\sigma$ is large enough, that $(\mathbb{L})$ admits a full basis of rational solutions. We proceed as in §8.4.1 to prove $\mathcal{A}^{[\sigma]} = \mathbb{L}^N$ and to identify a candidate rational vector solution $K^T \in \mathbb{L}(x)^{N+1}$. If $W^{[\sigma]}(z, 1)^T = 0$, then $K$ must be proportional to $(z, 1)^T$ in $\mathbb{L}((x))^{N+1}$, hence each series $z_i$ is equal to the rational series $K_i/K_{N+1}$. We therefore consider a parametrized candidate in the form $Q/P = \sum_{i=1}^N g_i(K_i/K_{N+1})$ for the parameter $a = g$ in $\mathbb{L}^N$.

* Else, $\rho^{[\sigma]} \le N-1$, and we proceed as in §8.4.2. The characteristic cone $\mathcal{A}^{[\sigma]}$ is given by the vanishing of all minors of rank $\rho^{[\sigma]} + 2$ in $W_a^{[\sigma]}$. (These are the minors of rank $\rho^{[\sigma]} + 1$ in the submatrix obtained by removing the last row and last column.) Extracting coefficients with regard to $x$ yields a linear system in $a$, so the cone $\mathcal{A}^{[\sigma]}$ is a vector space that can be parametrized in the form $a = gS$. Solving for a left kernel after this specialization delivers a basis of rows of the form $(\Pi_1, \dots, \Pi_{\rho^{[\sigma]}}, P, -Q)$, each parametrized by $g \in \mathbb{L}^v$. If $W^{[\sigma]}(z, 1)^T = 0$, then for each row

<footnote n="4">⁴We are indebted to Alin Bostan for asking us a question that led to the present section.</footnote>