# Quasi-Linear Guessing of Minimal Lexicographic Gröbner Bases of Ideals of C-Relations of Random Bi-Indexed Sequences

Jérémy Berthomieu, Romain Lebreton, Kevin Tran

► To cite this version:

Jérémy Berthomieu, Romain Lebreton, Kevin Tran. Quasi-Linear Guessing of Minimal Lexicographic Gröbner Bases of Ideals of C-Relations of Random Bi-Indexed Sequences. ISSAC '25: International Symposium on Symbolic and Algebraic Computation, Jul 2025, Guanajuato, Mexico. pp.258-266, 10.1145/3747199.3747569 . hal-04937888v2

HAL Id: hal-04937888

https://hal.sorbonne-universite.fr/hal-04937888v2

Submitted on 16 Jun 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
# Quasi-Linear Guessing of Minimal Lexicographic Gröbner Bases of Ideals of C-Relations of Random Bi-Indexed Sequences

Jérémy Berthomieu

jeremy.berthomieu@lip6.fr

Sorbonne Université, CNRS, LIP6

Paris, France

Romain Lebreton

Université de Montpellier, CNRS, LIRMM

Montpellier, France

romain.lebreton@lirmm.fr

Kevin Tran

Sorbonne Université, CNRS, LIP6

Paris, France

kevin.tran@lip6.fr

## ABSTRACT

Computing recurrence relations for sequences is a central problem in computer algebra, with applications in error-correcting codes, Gröbner basis computation, and sparse interpolation. While uni-indexed C-recursive sequences benefit from quasi-linear algorithms leveraging the half-gcd method, the extension to multi-indexed sequences remains computationally challenging. Existing methods for bi-indexed sequences achieve quadratic complexity at best, limiting their practical use.

This paper presents a quasi-linear algorithm for computing minimal lexiconographic Gröbner bases of the ideal of C-relations associated to bi-indexed sequences. Our approach extends the half-gcd algorithm in $K^N[y]$ by integrating a pseudo-Euclidean division. This approach shows how to leverage the bi-Hankel structure of the matrix, significantly improving the efficiency of computing C-relations and closing the complexity gap between the uni- and bi-indexed cases. Our algorithm is restricted to bi-indexed sequences whose associated bi-Hankel matrix has generic row rank profile.

## CCS CONCEPTS

* Computing methodologies → Algebraic algorithms.

## KEYWORDS

multi-indexed sequences, linear relation guessing, Hankel matrices, Gröbner bases, half-gcd algorithm, quasi-linear algorithm

## ACM Reference Format:

Jérémy Berthomieu, Romain Lebreton, and Kevin Tran. 2025. Quasi-Linear Guessing of Minimal Lexicographic Gröbner Bases of Ideals of C-Relations of Random Bi-Indexed Sequences. In *International Symposium on Symbolic and Algebraic Computation (ISSAC '25), July 28-August 1, 2025, Guanajuato, Mexico*. ACM, New York, NY, USA, 9 pages. https://doi.org/XXXXXXX. XXXXXXX

# 1 INTRODUCTION

*Context.* Guessing the minimal linear recurrence relation with constant coefficients (C-relation) of order *d* of a sequence $(u_i)_{i \in \mathbb{N}}$

is a fundamental problem in computer algebra and error correcting codes. It is for instance one of the latter steps of the Wiedemann algorithm [27] for computing the minimal polynomial of a matrix or solving a sparse linear system. The multi-indexed analogue, that is with a sequence $(u_{i_1,...,i_n})_{(i_1,...,i_n) \in \mathbb{N}^n}$ is at the root of *n*-dimensional cyclic codes and also the SPARSE-FGLM variant [12] of the FGLM algorithm [11] for Gröbner bases change of order.

Given the $D+1$ first terms $u_0, \dots, u_D$ of a uni-indexed sequence, the problem of computing the minimal C-relation can be modeled through a kernel computation of a Hankel matrix. It computes the correct relation as long as $D \ge 2d$. This Hankel structure leads to non-naive algorithms with complexity much better than $O(D^\omega)$, where $2 \le \omega < 3$ is the matrix multiplication exponent, relying on the extended Euclidean algorithm called on polynomials $x^{D+1}$ and $\sum_{i=0}^{D} u_i x^{D-i}$. The first instance of such a non-naive algorithm is due independently to Berlekamp [2] and Massey [17], both targeting an application to error correcting codes, and is now known as the BERLEKAMP-MASSEY algorithm. Thanks to quasi-linear algorithms for computing the extended Euclidean algorithm [9, 18], see also [13, Chap. 11], the complexity of computing such a minimal C-relation of order *d* is $\tilde{O}(D)$, as long as $D \ge 2d$.

*Related work.* The case of an *n*-indexed sequence, $n \ge 2$, $u = (u_{i_1,...,i_n})_{(i_1,...,i_n) \in \mathbb{N}^n}$ is more involved. The set of relations of $u$ forms an ideal, denoted $I(u)$, which is 0-dimensional whenever $u$ is C-recursive. Guessing consists in computing a representation of this ideal, that is a <-Gröbner basis for a given monomial order $<$ in the context of this paper. Denoting $S_u$ the <-staircase of $I(u)$, i.e. the monomials that are not <-leading monomials of $I(u)$, and $G_u$ the <-reduced Gröbner basis of $I(u)$, the complexity of the problem must depend on the number of given terms of $u$, and on $|G_u|$ and $|S_u|$, in order to encode the output in the monomial basis. The first algorithm to guess such a Gröbner basis is due to Sakata and extends the BERLEKAMP-MASSEY algorithm, leading the author to calling it the BERLEKAMP-MASSEY-SAKATA algorithm [21–23]. More recent algorithms were proposed based on linear algebra, i.e. computing the kernel of a multi-Hankel matrix [3, 4] or using a Gram-Schmidt process [19]. Another approach is based on multivariate polynomial arithmetic, especially division of polynomials such as [5, 6], or specifically for the bivariate case [14] using an approach similar to the uni-indexed case as they work on the polynomial $\sum_{j=0}^{D_y} (u_{i,j})_{i \in \mathbb{N}} y^{D_y-j} \in K^N[y]$. Finally, let us mention a bivariate Padé approximation method [20].

The complexity analysis of all these algorithms is not an easy task. Restricting ourselves to the case where the number of known terms of $u$ is minimal to ensure the correctness of the output allows us to express their complexities more easily. In the uni-indexed case, this would imply $D = \Theta(d)$, so that the complexity is $\tilde{O}(d)$.

The authors are supported by the joint ANR-FWF ANR-22-CE91-00075 EAGLES project and the ANR-21-CE39-0006 SANGRIA and ANR-22-PECY-0010 CRYPTANALYSE projects. We thank the referees for their valuable comments on the paper.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from permissions@acm.org.

ISSAC '25, July 28th - August 1st, 2025, Guanajuato, Mexico

© 2025 Copyright held by the owner/author(s). Publication rights licensed to ACM.

ACM ISBN 978-1-4503-XXXXX-X/2018/06...$15.00

https://doi.org/XXXXXXX.XXXXXXX
In [23], the complexity of the BERLEKAMP-MASSEY-SAKATA algorithm is $O(|S_u|^2 \cdot |G_u|)$, though the output need not be a reduced Gröbner basis. The complexity of the algorithm of [3, 4] is $O(|S_u|^\omega + |S_u|^2 \cdot |G_u|)$ and the output is reduced. The algorithm of [19] has complexity $O(|S_u|^2 \cdot (|S_u| + |G_u|))$, where $G_u$ is a border basis, and thus has larger size than $G_u$, while the algorithm [5, 6] has a similar complexity $O(|S_u|^2 \cdot (|S_u| + |G_u|))$. Furthermore, they all need the sequence terms $u_{i_1, \dots, i_n}$ where $x_i^{i_1} \cdots x_n^{i_n}$ is in the Minkowski sum of $S_u$ with itself, denoted $2S_u$. If we simplify further to the bi-indexed case and we denote $d_x$ (resp. $d_y$) the maximal degree in $x$ (resp. $y$) of $G_u$, these complexity upper bounds become at least $O(\max(d_x, d_y)^2 |G_u|)$ using the fact that $|S_u| \ge d_x + d_y - 1$. Now, on the one hand, all the monomials $x^i y^j$ for $0 \le i < d_x$ and $0 \le j < d_y$ are in $2S_u$ and, on the other hand, all monomials in $2S_u$ have degree in $x$ (resp. $y$) at most $2d_x - 2$ (resp. $2d_y - 2$). Hence, all these algorithms need exactly $\Theta(d_x d_y)$ terms on input. Finally, using also $\Theta(d_x d_y)$ terms of $u$, the algorithm of [14] computes a Gröbner basis of $I(u)$ in $\tilde{O}(d_x^{\omega+1} d_y)$ operations, while [20] requires $\tilde{O}(\min(d_x, d_y)^\omega d_x d_y)$ operations.

**Contribution.** The main contribution of this paper is GUESSING-BIVAR, an algorithm that takes as input the $(D_x + 1)(D_y + 1)$ sequence terms $u_{i,j}$ for $0 \le i \le D_x$ and $0 \le j \le D_y$ and returns a minimal lexicographic Gröbner basis of $I(u)$, with support in $\{x^i y^j | 0 \le i \le d_x, 0 \le j \le d_y\}$ for $u = (u_{i,j})_{(i,j) \in \mathbb{N}^2}$ using $\tilde{O}(D_x D_y + d_x d_y |G_u|)$ operations. This algorithm works under the assumption that the multi-Hankel matrix $(u_{i+k,j+\ell})_{x^i y^j, x^k y^\ell} \in S_u$ has an LU decomposition without pivoting. This condition is experimentally always satisfied whenever the terms $(u_{i,j})_{x^i y^j} \in S_u$ are picked at random. As a consequence, this closes the complexity gap between the uni-indexed case and the bi-indexed one.

**Motivations.** The study of guessing algorithm is motivated by polynomial system solving application. In future works, we would like to adapt the SPARSE-FGLM algorithm for specific non shape position cases, such as 2-thick ones [1]. Furthermore, guessing is related to structured system solving and to computing the bivariate resultant, see [25, 26].

**Organization of the paper.** In §2, we recall the polynomial representation of C-relations, and also how to relate their guessing to linear algebra and univariate gcd computation. In §3, we extend this viewpoint to bi-indexed sequences under the aforementioned assumption on the associated multi-Hankel matrix. In §4, we design a half-gcd-like algorithm on bivariate polynomials and how it can be used as a subroutine of GUESSINGBIVAR for guessing. Finally, our benchmarks in §5 confirm the efficiency of our algorithm.

# 2 PRELIMINARIES

In this section, we recall all basic definitions and results on C-recursive uni- and bi-indexed sequences, polynomials and Gröbner bases. We consider $\mathbb{N}$ as the set of all natural numbers including 0, also consider that $\deg(0) = -\infty$. If there is no ambiguity on the number of indices we denote the sequence as $u$.

## 2.1 Uni-indexed sequences

For uni-indexed sequences, C-recursive sequences are the ones satisfying linear recurrences with constant coefficients. A sequence $(u_i)_{i \in \mathbb{N}}$ is C-recursive if there exist $g_0, \dots, g_{d-1} \in \mathbb{K}$ such that

for $i \in \mathbb{N}$, $u_{i+d} = g_{d-1}u_{i+d-1} + \cdots + g_0u_i$. Such a combination is called C-relation and can be represented as a polynomial $g = x^d - \sum_{i=0}^{d-1} g_i x^i \in \mathbb{K}[x]$. Computing a C-relation can be reduced to a linear system solving problem.

The Hankel matrix of size $d$ associated to the sequence $u = (u_i)_{i \in \mathbb{N}}$ is $\mathcal{H} = (u_{i+k})_{0 \le i, k < d} \in \mathbb{K}^{d \times d}$. Moreover, one can compute the C-relation $g = x^d - \sum_{i=0}^{d-1} g_i x^i$ by solving the linear system $[g_0 \ g_1 \ \cdots \ g_{d-1}]^\mathcal{H} = [u_d \ \cdots \ u_{2d-1}]$.

Another approach is done using generating series $S = \sum_{i \in \mathbb{N}} u_i x^i \in \mathbb{K}[x]$, as the generating series of a C-recursive sequence admits a finite representation. Indeed, for such series $S$, there exists $p, q \in \mathbb{K}[x]$ such that $qS = p$ with $\deg(p) < d$ and $q = \hat{g} := x^d g(1/x) \in \mathbb{K}[x]$ is the reversal of $g$.

From the degree constraint on $p$ and $q$, one can recover $p$ and $q$ from the relation $qS = p$ mod $x^{D+1}$ with $D \ge 2d$. This modular equation can be rewritten as a Bézout's identity $qS + rx^{D+1} = p$ with $r \in \mathbb{K}[x]$ and computing $q = \hat{g}$ comes down to computing a Truncated Extended Euclidean algorithm. A fast computation of this relation can be done through a call to the half-gcd algorithm [15, 18, 24].

**LEMMA 2.1.1.** Let $a, b \in \mathbb{K}[x]$ with $\deg(a) = D$ and $\deg(b) = d$ such that $D \ge d$. Computing $q, r \in \mathbb{K}[x]$ satisfying $a =qb+r$ with $\deg(r) < d$ can be done in $\tilde{O}(D)$ operations in $\mathbb{K}$.

The transpose of this operation, called the extension, is computed in the same complexity by the Tellegen principle [8]. This operation corresponds to the extension of C-recursive sequences $(u_i)_{i \in \mathbb{N}}$ by the C-relation $g = x^d - \sum_{i=0}^{d-1} g_i x^i$, for $S = \sum_{i=0}^{d-1} u_i x^i$, it computes $\tilde{S} = \sum_{i=0}^{D} u_i x^i$ using $u_{i+d} = g_{d-1}u_{i+d-1} + \cdots + g_0u_i$ for $i \ge 0$.

Recall that $\tilde{O}(\cdot)$ means that polylogarithmic factors are omitted.

**THEOREM 2.1.2.** Computing the C-relation $g \in \mathbb{K}[x]$ on $(u_i)_{i \in \mathbb{N}}$ of degree $d$, knowing the $D+1$ initial terms of $(u_i)_{i \in \mathbb{N}}$ with $D \ge 2d$, can be done in $\tilde{O}(D)$ operations in $\mathbb{K}$.

The half-gcd algorithm can also be derived to a Hankel system solving of size $d-1$ (see [9]) and can be done in $\tilde{O}(d)$ operations.

## 2.2 Bivariate polynomial rings

For bi-indexed sequences, we use bivariate polynomials to represent the C-relations. For a polynomial $f \in \mathbb{K}[x, y]$ and $(i, j) \in \mathbb{N}^2$, we note $f_{i,j}$ the coefficient of $f$ associated to the monomial $x^i y^j$, the support of $f$ is the monomial set supp($f$) = $\{x^i y^j | f_{i,j} \ne 0\}$.

We define the box monomial set of parameter $(d_x, d_y) \in \mathbb{N}^2$ as $(d_x, d_y)$-box := $\{x^i y^j | 0 \le i \le d_x \text{ and } 0 \le j \le d_y\}$. We note by $\mathbb{K}[x, y]_{\le (d_x, d_y)}$ the set of polynomials with support in $(d_x, d_y)$-box.

For polynomials $f, g \in \mathbb{K}[x, y]_{\le (d_x, d_y)}$, the addition of $f+g$ can be computed using $O(d_x d_y)$ operations in $\mathbb{K}$ and the multiplication $fg$ can be computed using $\tilde{O}(d_x d_y)$ operations.

For bivariate polynomials, we have to define a total order on the monomial set. In our study, we are only interested in the lexicographic order, we refer to [10] for more general consideration. We note $<$ the lexicographic order on $\mathbb{K}[x, y]$ with $x < y$ and such that $x^i y^j < x^k y^\ell$ if $j < \ell$, $j = \ell$ and $i < k$. For a nonzero polynomial $f \in \mathbb{K}[x, y]$, the leading monomial $\text{lm}(f)$ of $f$ w.r.t. $<$ is the maximum monomial of $f$ ordered by $<$. The leading coefficient $\text{lc}(f)$ of $f$ w.r.t. $<$ is the coefficient associated to $\text{lm}(f)$. The leading term $\text{lt}(f)$ of $f$ w.r.t. $<$ satisfies $\text{lt}(f) = \text{lc}(f) \text{lm}(f)$.
An ideal of $K[x, y]$ can be generated by a finite set of polynomials. Gröbner bases are particular sets of generators with interesting computational properties. For an ideal $I \subseteq K[x, y]$, a Gröbner basis $\mathcal{G}$ of $I$ for the lexicographic order is a finite generating set of $I$ such that $\langle \mathrm{lm}(\mathcal{G}) \rangle = \langle \mathrm{lm}(I) \rangle$, i.e. it spans $\mathrm{lm}(I)$ as a monomial set. A minimal Gröbner basis $\mathcal{G}$ is a Gröbner basis such that no $\mathrm{lm}(g) \in \mathrm{lm}(\mathcal{G})$ is divisible by an element in $\mathrm{lm}(\mathcal{G}\setminus\{g\})$ and for all $g \in \mathcal{G}$, $\mathrm{lc}(g) = 1$. The (unique) reduced Gröbner basis $\mathcal{G}$ is a minimal Gröbner basis such that for all $g \in \mathcal{G}$, the monomials $m \in \mathrm{supp}(g)$ are not divisible by any $\mathrm{lm}(\mathcal{G}\setminus\{g\})$.

The staircase $S$ associated to an ideal $I$ is $S := \{x^i y^j \mid x^i y^j \notin \mathrm{lm}(I)\}$. It forms a $K$-vector space basis of the quotient ring $K[x, y]/I$. The polynomial division with remainder by a Gröbner basis (defined in [10, Chapter 2.7]) gives a unique polynomial $r$ with support in the staircase $S$. For $f \in K[x, y]$, we denote by $r = f \bmod \mathcal{G} \in K[x, y]$ with $\mathrm{supp}(r) \subset S$ the unique remainder of $f$ by a Gröbner basis $\mathcal{G}$ [10, Chapter 2.7]. Let $I$ be an ideal of $K[x, y]$ and let $f \in K[x, y]$, the colon ideal of $I$ by $f$ is $I : \langle f \rangle = \{g \in K[x, y] \mid gf \in I\}$ (see [10, Chapter 4.4] for a more general description).

## 2.3 Bi-indexed C-recursive sequences

The set $K^{\mathbb{N}^2}$ corresponds to the set of bi-indexed sequence with terms in $K$. We denote the zero sequence by $0 = (0)_{(i,j)\in\mathbb{N}^2}$. For C-recursive sequences, we allow two types of operations: index shifts and scalar multiplications on sequences. These operations can be described by the action $\star$ of $K[x, y]$ in $K^{\mathbb{N}^2}$ such that $x^k y^l \star u = (u_{i+k,j+l})_{(i,j)\in\mathbb{N}^2}$ for $k, l \in \mathbb{N}$, and extended by linearity to $K[x, y]$.

For a sequence $u$, a C-relation $g$ on $u$ is a polynomial $g \in K[x, y]$ that satisfies $g \star u = 0$. We note by $I(u) = \{g \in K[x, y] \mid g \star u = 0\}$ the ideal of relations of $u$. A sequence $u$ is C-recursive if the ideal of relations $I(u)$ is 0-dimensional i.e. $\dim_K(K[x, y]/I(u)) < \infty$.

Furthermore, we denote by $\mathcal{G}_u$ the reduced $\langle$-Gröbner basis of $I(u)$ and by $S_u$ the $\langle$-staircase of $I(u)$. We also let $S_{u, <m} = \{x^i y^j \in S_u \mid x^l y^j < m\}$. We note the exponents set of $S_u$ by $\mathcal{E}_u = \{(i, j) \in \mathbb{N}^2 \mid x^i y^j \in S_u\}$ and $\mathcal{E}_{u, <(k,l)} = \{(i, j) \in \mathcal{E}_u \mid x^i y^j < x^k y^l\}$.

For $u$ a C-recursive sequence and $\mathcal{G}$ a Gröbner basis of $I(u)$, any term of $u$ can be computed from the relations in $\mathcal{G}$ and the initial terms in $S_u$ [22]. A C-recursive sequence is uniquely determined by the terms associated to $\mathcal{E}_u$, as the other terms are linear combinations thereof.

**LEMMA 2.3.1 ([23, §2]).** Fix $I(u)$ and $\mathcal{G}$ a Gröbner basis of $I(u)$ w.r.t. the order $<$. Then $I(u) \subseteq I(v)$ iff for all $(k, l) \in \mathrm{lm}(I(u))$, we have $v_{k,l} = \sum_{(i,j) \in \mathcal{E}_u} c_{i,j} v_{i,j}$ with $x^k y^l \bmod \mathcal{G} = \sum_{(i,j) \in \mathcal{E}_u} c_{i,j} x^i y^j$.

We define the $K$-linear subspace $L_u := \{h \star u \mid h \in K[x, y]\} \subset K^{\mathbb{N}^2}$ and consider the linear application $\phi(h) = h \star u$ from $K[x, y]$ to $L_u$. By construction, $\phi$ is surjective. As ker $\phi = I(u)$, we can define the isomorphism $\bar{\phi} : K[x, y]/I(u) \to L_u$ from $\phi$. We define $\mathcal{F} = \{e_{i,j}\}_{(i,j)\in\mathcal{E}_u} \subset K^{\mathbb{N}^2}$ with $e_{i,j}$ defined for $(k, l) \in \mathcal{E}_u$ such that $(e_{i,j})_{k,l} = 0$ if $(i, j) \neq (k, l)$ and $(e_{i,j})_{i,j} = 1$ and outside $\mathcal{E}_u$ we extend the terms of $e_{i,j}$ in $\mathcal{E}_u$ with the relations in $I(u)$.

We note $\mathcal{H}_{S_u}$ the matrix associated to $\bar{\phi}$, with the basis $S_u$ for $K[x, y]/I(u)$ and $\mathcal{F}$ for $L_u$ both ordered w.r.t. $<$. The application $\bar{\phi}$ is an isomorphism so the matrix $\mathcal{H}_{S_u}$ is invertible.

**THEOREM 2.3.2.** Let $u$ and $v$ be two C-recursive sequences. The following statements are equivalent:

(a) $\exists! h \in K[x, y]$ with support in $S_u$ such that $v = h \star u$;

(b) $\exists h \in K[x, y]$ such that $I(v) = I(u) : \langle h \rangle$;

(c) $I(u) \subset I(v)$.

*PROOF.* For (a) $\Rightarrow$ (b), we have $g \in I(v) \Leftrightarrow 0 = g \star v = (gh) \star u \Leftrightarrow g \in I(u) : \langle h \rangle$. For (b) $\Rightarrow$ (c), it is direct by definition.

For (c) $\Rightarrow$ (a), let's first prove that $I(u) \subset I(v)$ implies $v \in \mathrm{span}_K(F)$. Let $w = v - \sum_{(i,j) \in \mathcal{E}_u} v_{i,j} e_{i,j}$. For $(k, l) \in \mathcal{E}_u$, we have by construction $w_{k,l} = 0$. From Lm. 2.3.1, we have $I(u) \subset I(e_{i,j})$ for any $(i, j) \in \mathcal{E}_u$. Since $f \in I(u)$ is in $I(v)$ and all $I(e_{i,j})$, we deduce that $I(u) \subset I(w)$. Hence, $w = 0$ and $v = \sum_{(i,j) \in \mathcal{E}_u} v_{i,j} e_{i,j}$.

We now prove that the family $\mathcal{F}$ is a basis of $L_u$. By construction, $\mathcal{F}$ is linearly independent. Let $h \star u \in L_u$, since for $f \in I(u)$, $(fh) \star u = f \star (h \star u) = 0$, we have $I(u) \subset I(h \star u)$. As seen before, this implies $h \star u \in \mathrm{span}_K(F)$, so $L_u \subset \mathrm{span}_K(F)$. Since $\dim_K(L_u) = \dim_K(K[x, y]/I(u)) = |S_u|$, we conclude that $\mathcal{F}$ is a basis of $L_u$. Therefore, we can write $v = \sum_{(i,j) \in \mathcal{E}_u} v_{i,j} e_{i,j}$ so $v \in L_u$.

For uniqueness, let $h' \in K[x, y]$ s.t. $v = h' \star u$ and $\mathrm{supp}(h') \subset S_u$. We get $(h - h') \star u = 0$ so $h - h' \in I(u)$ and $h - h' \mathrm{rem} \mathcal{G}_u = 0$. Since $\mathrm{supp}(h)$, $\mathrm{supp}(h') \subset S_u$ by the linearity of the reduction $h = h \mathrm{rem} \mathcal{G}_u = h' \mathrm{rem} \mathcal{G}_u = h'$. $\square$

# 3 BI-INDEXED SEQUENCES

In this section, we restrict ourselves to C-recursive bi-indexed sequences $u = (u_{i,j})_{i,j \in \mathbb{N}}$. We denote by $d_x, d_y \in \mathbb{N}$ the exponents satisfying $x^{d_x}, y^{d_y} \in \mathrm{lm}(\mathcal{G}_u)$.

## 3.1 Hankel matrix and LU decomposition

For a bi-indexed sequence $u$ and $j \in \mathbb{N}$, we note the sub-sequences $u_{*,j} = (u_{i,j})_{i \in \mathbb{N}} \in K^{\mathbb{N}}$. Sub-sequences do not necessarily contain enough information to recover the ideal $I(u) \cap K[x]$. For example, let $u = ((-1)^{ij})_{i,j \in \mathbb{N}}$, then $I(u_{*,j}) = \langle x - (-1)^j \rangle$, but $I(u) \cap K[x] = \langle x^2 - 1 \rangle$. To overcome the problem posed by this example, we make the following assumption on the sequence $u$.

**Assumption A.** The matrix $\mathcal{H}_{S_u}$ defined in §2.3 for the bi-indexed sequence $u$ admits an LU decomposition.

For a matrix $M \in K^{n \times n}$, the principal $r \times r$ submatrix $M_r \in K^{r \times r}$ is the matrix built from the first $r$ rows and columns of $M$. Recall that an invertible matrix $M \in K^{n \times n}$ admits an LU decomposition iff for $1 \le r \le n$, the submatrix $M_r$ is invertible.

Consider $\mathcal{H}_{S_u} = LU$ with $L$ a lower triangular matrix with ones on the diagonal and $U$ an upper triangular matrix. We note the rows of $L^{-1} = [\ell_m]_{m \in S_u}$ with $\ell_m \in K^{1 \times |S_u|}$. We note $p_m \in K[x, y]$ the polynomial representing $\ell_m$ in the basis $S_u$. The matrix $L^{-1}$ is lower triangular with ones on its diagonal, so $\mathrm{lt}(p_m) = m$ for all $m \in S_u$. For $y^j \in S_u$, we denote by $u^{(j)}$ the sequence $p_{y^j} \star u$.

**LEMMA 3.1.1.** For $0 \le j < d_y$, we have $u_{*,k}^{(j)} = 0$ for $0 \le k < j$.

*PROOF.* Let $i \in \mathbb{N}$, $k < j$ and consider the term $u_{i,k}^{(j)}$ of $u^{(j)}$. By construction, the row of $U$ indexed by $y^j$ contains terms of $u^{(j)}$ and in particular $(u^{(j)})_{r,s} = 0$ for $(r, s) \in \mathcal{E}_{u, <(0,j)}$. Now, since $u_{i,k}^{(j)} = (x^i y^k \star u^{(j)})_{0,0} = ((x^i y^k \mathrm{rem} \mathcal{G}_u) \star u^{(j)})_{0,0}$, we express $u_{i,k}^{(j)}$ as a linear combination of $u_{r,s}^{(j)} = 0$ for $(r, s) \in \mathcal{E}_{u, <(0,j)}$. $\square$

**LEMMA 3.1.2.** Let $j \in \mathbb{N}$ and $t \in K[x, y]$ such that $(t \star u)_{r,s} = 0$ for $(r, s) \in \mathcal{E}_{u, <(0,j)}$. If $\mathrm{deg}_y(t) < j$, then $t \in I(u)$. Otherwise, if $\mathrm{lt}(t) = y^j$ for $0 \le j < d_y$, then $p_{y^j} = t \mathrm{rem} \mathcal{G}_u$.
PROOF. Let $\mathcal{H}_{S_{u,<yj}}$ be the principal submatrix of $\mathcal{H}_{S_u}$ with rows indexed by $S_{u,<yj}$ and columns by $\mathcal{E}_{u,<(0,j)}$. If $\deg_y(t) < j$, we can represent the polynomial $\tilde{t} := t \text{ rem } \mathcal{G}_u$ by a vector $\ell$ in the basis $S_{u,<yj}$. Since $(\tilde{t} \star u)_{r,s} = 0$ for $(r,s) \in \mathcal{E}_{u,<(0,j)}$, $\ell$ satisfies $\ell \mathcal{H}_{S_{u,<yj}} = 0$. As $u$ satisfies Asm. A, $\mathcal{H}_{S_{u,<yj}}$ is invertible so $\ell = 0$, and $t \in I(u)$. Now, if $\text{lt}(t) = y^j$ for $0 \le j < d_y$, then $\tilde{t} := t - p_{yj}$ satisfies the hypotheses and $\deg_y(\tilde{t}) < j$, so $\tilde{t} \in I(u)$ and $p_{yj} = p_{yj}$ \text{rem } \mathcal{G}_u = t \text{ rem } \mathcal{G}_u$. $\square$

THEOREM 3.1.3. The family $\mathcal{P} = \{x^i p_{yj} | x^i y^j \in S_u\}$ is a basis of $\mathbb{K}[x, y]/I(u)$ as a $\mathbb{K}$-vector space.

PROOF. Let $x^i y^j \in S_u$ and consider $g = x^i p_{yj}$ \text{rem } \mathcal{G}_u$. Since $\text{lt}(x^i p_{yj}) = x^i y^j \in S_u$, we have $g = x^i y^j + \sum_{(r,s) \in \mathcal{E}_{u,<(i,j)}} c_{r,s} x^r y^s$. Hence, the change-of-basis matrix between the bases $S_u$ and $\mathcal{P}$ is lower triangular with ones on its diagonal. $\square$

We define the matrix $\mathcal{H}_{\mathcal{P}}$ representing the application $\bar{\phi}$ (see §2.3) with row basis $\mathcal{P}$ and column basis $\mathcal{F}$ defined in §2.3.

LEMMA 3.1.4. The matrix $\mathcal{H}_{\mathcal{P}}$ is block upper triangular and its diagonal blocks are invertible, i.e.

$$ \mathcal{H}_{\mathcal{P}} = \begin{bmatrix} \mathcal{H}_0 & \mathcal{H}_1 & \cdots & \mathcal{H}_{d_y-1} \\ 0 & \mathcal{H}_1^{(1)} & \cdots & \mathcal{H}_{d_y-1}^{(1)} \\ \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathcal{H}_{d_y-1}^{(d_y-1)} \end{bmatrix} $$

PROOF. For $0 \le j < d_y$, the $j$th row block starts with $j$ zero matrices since $\mathbf{u}_{*,k}^{(j)} = 0$ for $0 \le k < j$ from Lm. 3.1.1 thus the matrix $\mathcal{H}_{\mathcal{P}}$ is block upper triangular matrix. The matrix $\mathcal{H}_{\mathcal{P}}$ is invertible since the linear application $\bar{\phi}$ is an isomorphism, hence the block diagonal matrices $\mathcal{H}_j^{(j)}$ are invertible. $\square$

THEOREM 3.1.5. For $0 \le j < d_y$ and $d_j \in \mathbb{N}$, there exists $g \in I(u)$ such that $\text{lm}(g) = x^{d_j} y^j$ if and only if there exists $f_j \in I(u^{(j)}) \cap \mathbb{K}[x]$ such that $\text{lm}(f_j) = x^{d_j}$.

PROOF. Let $f_j \in I(u^{(j)}) \cap \mathbb{K}[x]$ with $\text{lm}(f_j) = x^{d_j}$, by definition $0 = f_j \star u^{(j)} = (f_j p_{yj}) \star u$ so $g = f_j p_{yj} \in I(u)$ and $\text{lm}(g) = x^{d_j} y^j$. $\square$

Let $g \in I(u)$ with $\text{lm}(g) = x^{d_j} y^j \notin S_u$ by definition. Consider the sequence $(x^{d_j} p_{yj}) \star u = x^{d_j} \star u^{(j)}$ from Lm. 3.1.1 we have for $k < j$, $(x^{d_j} \star u^{(j)})_{*,k} = 0$. From Lm. 3.1.4, the matrix $\mathcal{H}_j^{(j)}$ is invertible, so there exists a polynomial $f \in \mathbb{K}[x]$ with $\text{supp}(fy^j) \subset S_u$ satisfying $(f \star u^{(j)})_{r,j} = (x^{d_j} \star u^{(j)})_{r,j}$ for $(r,j) \in \mathcal{E}_u$. By construction, the polynomial $t = (x^{d_j} - f)p_{yj}$ is such that $(t \star u)_{r,s} = 0$ for $(r,s) \in \mathcal{E}_{u,<(0,j+1)}$ and $\deg_y(t) = j$. So by Lm. 3.1.2, $t \in I(u)$ and $f_j = x^{d_j} - f \in I(u^{(j)})$. Note that $x^{d_j} y^j \notin S_u$ and $\text{supp}(fy^j) \subset S_u$, so $\text{lm}((x^{d_j} - f)p_{yj}) = \text{lm}(x^{d_j} p_{yj}) = x^{d_j} y^j$ and $\text{lm}(f_j) = x^{d_j}$. $\square$

THEOREM 3.1.6. For $0 \le j < d_y$, the sequence $\mathbf{u}_{*,j}^{(j)} \in \mathbb{K}^N$ is such that $I(\mathbf{u}_{*,j}^{(j)}) = \bigcap_{k \ge j} I(\mathbf{u}_{*,k}^{(j)})$ i.e. $\mathbf{u}_{*,k}^{(j)} \in L_{\mathbf{u}_{*,j}^{(j)}}$ for $k \ge j$.

PROOF. The inclusion $\bigcap_{k \ge j} I(\mathbf{u}_{*,k}^{(j)}) \subset I(\mathbf{u}_{*,j}^{(j)})$ is direct. For the reverse inclusion, let $g \in I(\mathbf{u}_{*,j}^{(j)})$. The polynomial $t = gp_{yj}$ is such

that $(t \star u)_{r,s} = 0$ for $(r,s) \in \mathcal{E}_{u,<(0,j+1)}$, since it is zero for $s < j$ by Lm. 3.1.1, and for $s = j$ by definition of $g$. Since $\text{deg}_y(t) \le j$ (as $\text{lt}(p_{yj}) = y^j$), we have $t \in I(u)$ using Lm. 3.1.2, and thus $g \star u^{(j)} = t \star u = 0$ so $g \in \bigcap_{k \ge j} I(\mathbf{u}_{*,k}^{(j)})$. $\square$

From this property, we can find a relation between the sequences $(\mathbf{u}^{(j)})_{0 \le j \le d_y}$ with $\mathbf{u}^{(d_y)} = 0$.

THEOREM 3.1.7. For $0 \le j < d_y$, we have $\mathbf{u}_{*,j+1}^{(j+1)} \in L_{\mathbf{u}_{*,j}^{(j)}}$ and there exists $(\bar{a}_j, \bar{b}_j) \in \mathbb{K}[x]^2$ with $\text{supp}(\bar{a}_j) \subset \mathcal{S}_{\mathbf{u}_{*,j-1}^{(j-1)}}$ and $\text{supp}(\bar{b}_j) \subset \mathcal{S}_{\mathbf{u}_{*,j}^{(j)}}$ satisfying $\mathbf{u}^{(j+1)} = \bar{a}_j \star \mathbf{u}^{(j-1)} + (y - \bar{b}_j) \star \mathbf{u}^{(j)}$ if $j \ne 0$ and $\mathbf{u}^{(j)} = (y - \bar{b}_0) \star \mathbf{u}$ if $j = 0$. $\square$

PROOF. We prove the statement by induction on $j$. For $j = 0$, we have by definition $\mathbf{u}^{(1)} = p_y \star u$ with $\text{lt}(p_y) = y$ and $\text{supp}(p_y) \subset S_u$ so $p_y = y - \bar{b}_0$ with $\text{supp}(\bar{b}_0) \subset S_{u,0}$. From this relation, we deduce that $\mathbf{u}_{*,1}^{(1)} = \mathbf{u}_{*,0} - \bar{b}_0 \star \mathbf{u}_{*,1}$ and by Thm. 3.1.6 it results that $\mathbf{u}_{*,1}^{(1)} \in L_{u,0}$. For $1 \le j < d_y - 1$, we suppose that the statement is true at step $j-1$ and prove that it holds at step $j$. Consider $(\bar{a}_j, \bar{b}_j) \in \mathbb{K}[x]^2$ with $\text{supp}(\bar{a}_j) \subset \mathcal{S}_{\mathbf{u}_{*,j-1}^{(j-1)}}$ and $\text{supp}(\bar{b}_j) \subset \mathcal{S}_{\mathbf{u}_{*,j}^{(j)}}$ satisfying $\bar{a}_j \star \mathbf{u}_{*,j-1}^{(j-1)} = -\mathbf{u}_{*,j}^{(j)}$ and $\bar{a}_j \star \mathbf{u}_{*,j}^{(j-1)} + \mathbf{u}_{*,j+1}^{(j)} = \bar{b}_j \star \mathbf{u}_{*,j}^{(j)}$. There exists $\bar{a}_j$ satisfying the first equality by the induction hypothesis $\mathbf{u}_{*,j}^{(j)} \in L_{\mathbf{u}_{*,j-1}^{(j-1)}}$. For $\bar{b}_j$, by Thm. 3.1.6 we have $\bar{a}_j \star \mathbf{u}_{*,j}^{(j-1)} \in L_{\bar{a}_j \star \mathbf{u}_{*,j-1}^{(j-1)}} = L_{\mathbf{u}_{*,j}^{(j)}}$ hence from Thm. 2.3.2 we can find $\bar{b}_j$ satisfying the conditions. Let $w = \bar{a}_j \star \mathbf{u}^{(j-1)} + (y - \bar{b}_j) \star \mathbf{u}^{(j)}$. By construction of $w$, we have $\mathbf{w}_{*,k} = 0$ for $k < j+1$ and $w = t \star u$ with $t = (\bar{a}_j p_{yj-1} + (y - \bar{b}_j) p_{yj})$. If $j \ne d_y - 1$ then by Lm. 3.1.2 since $\text{lt}(t) = y^{j+1}$ we deduce that $p_{yj+1} = t \text{ rem } \mathcal{G}_u$ and $\mathbf{u}^{(j+1)} = w$, otherwise if $j = d_y - 1$ then $\mathbf{w}_{*,k} = 0$ for $k \le d_y$ so $w = 0 = \mathbf{u}^{(d_y)}$. Finally, the relation $\mathbf{u}_{*,j+1}^{(j+1)} = \bar{a}_j \star \mathbf{u}_{*,j+1}^{(j-1)} + \mathbf{u}_{*,j}^{(j)} - \bar{b}_j \star \mathbf{u}_{*,j+1}^{(j)}$ gives $\mathbf{u}_{*,j+1}^{(j+1)} \in L_{\mathbf{u}_{*,j}^{(j)}}$ with the same arguments used to prove the existence of $\bar{b}_j$. $\square$

LEMMA 3.1.8. For $0 \le j < d_y$, if we define $[\begin{matrix} \bar{s}_j & \bar{t}_j \\ \bar{s}_{j+1} & \bar{t}_{j+1} \end{matrix}] = \bar{Q}_j \cdots \bar{Q}_0$ with $\bar{Q}_k = [\begin{matrix} 0 & 1 \\ \bar{a}_k & \bar{b}_k \end{matrix}]$ then $p_{yj} = \bar{t}_j \text{ rem } \mathcal{G}_u$ and $\bar{t}_{dy} \in I(u)$. $\square$

PROOF. For $0 \le j < d_y$, we have $[\begin{matrix} \mathbf{u}^{(j)} \\ \mathbf{u}^{(j+1)} \end{matrix}] = \bar{Q}_j \cdots \bar{Q}_1 [\begin{matrix} \mathbf{u}^{(0)} \\ \mathbf{u}^{(1)} \end{matrix}]$. If we note $R_j = \bar{Q}_j \cdots \bar{Q}_1 = [\begin{matrix} \alpha_j & \beta_j \\ \gamma_j & \delta_j \end{matrix}]$ then $[\begin{matrix} \bar{s}_j & \bar{t}_j \\ \bar{s}_{j+1} & \bar{t}_{j+1} \end{matrix}] = R_j [\begin{matrix} 0 & 1 \\ \bar{a}_0 & (y-\bar{b}_0) \end{matrix}]$ and $\bar{t}_j = \alpha_j + (y - \bar{b}_0)\beta_j$. From Thm. 3.1.7, we have $\mathbf{u}^{(1)} = (y - \bar{b}_0) \star \mathbf{u}^{(0)}$ so $\bar{t}_j \star \mathbf{u} = \alpha_j \star \mathbf{u} + \beta_j \star \mathbf{u}^{(1)} = \mathbf{u}^{(j)}$. By the same reasoning, we obtain $\bar{t}_{j+1} \star \mathbf{u} = \mathbf{u}^{(j+1)}$. Therefore, $(p_{yj} - \bar{t}_j) \star \mathbf{u} = 0$ so $(p_{yj} - \bar{t}_j \text{ rem } \mathcal{G}_u) = 0$ and by linearity of the reduction we get $p_{yj} = p_{yj} \text{ rem } \mathcal{G}_u = \bar{t}_j \text{ rem } \mathcal{G}_u$. For $j = d_y$, we have $\mathbf{u}^{(d_y)} = 0 = \bar{t}_{dy} \star \mathbf{u}$ so $\bar{t}_{dy} \in I(u)$. $\square$

## 3.2 Pseudo-Euclidean division

In this subsection, following [14, Sec. 6], we work with polynomials in $\mathbb{K}^N[y]$. We define an action $\cdot$ of $\mathbb{K}[x, y]$ on $\mathbb{K}^N[y]$ that mimics the action $\star$ on sequences. If $r = \sum_{j=0}^D r_j y^j \in \mathbb{K}^N[y]$, $g \in \mathbb{K}[x]$ and
$d \in \mathbb{N}$, we set $(gy^d) \cdot r = \sum_{j=0}^{D}(g \star r_j)y^{j+d}$, then extend linearly
· to $\mathbb{K}[x, y]$. As in the uni-indexed case, our goal is to reduce the
guessing problem to the computation of successive remainders. For
that we define a pseudo-Euclidean division in $\mathbb{K}^N[y]$.

THEOREM 3.2.1. Let $f = \sum_{j=0}^{d} f_j y^j$ and $g = \sum_{j=0}^{d-1} g_j y^j$ be two
polynomials in $\mathbb{K}^N[y]$ of respective degree $d$ and $d-1$ with $d \ge 1$.

If (i) $g_{d-1} \in L_{fd}$, (ii) $f_{d-1} \in L_{fd}$, (iii) $g_{d-2} \in L_{gd-1}$ with $fd$ C-
recursive then $\exists! (a, b) \in \mathbb{K}[x]^2$ with $\text{supp}(a) \subset S_{fd}$ and $\text{supp}(b) \subset
S_{gd-1}$ satisfying $a \cdot f = (-y+b) \cdot g + r$ with $\text{deg}_y(r) < \text{deg}_y(g)$.

When the conditions (i), (ii), (iii) of the previous theorem are
satisfied, we say that the pseudo-Euclidean division of $f$ by $g$ is
well-defined and that its result is $(a, b, r)$.

PROOF. Consider such polynomials $f, g \in \mathbb{K}^N[y]$, since $g_{d-1} \in L_{fd}$ from Thm. 2.3.2 there exists a unique polynomial $a \in \mathbb{K}[x]$ with $\text{supp}(a) \subset S_{fd}$ such that $g_{d-1} = -a \star f_d$. We can construct $\tilde{g} = a \cdot f + y \cdot g = \sum_{j=1}^{d-1} (a \star f_j + g_{d-1}) y^j + a \star f_0$. If $\text{deg}_y(\tilde{g}) < d-1$, then $\tilde{g} = r$ and the pair $(a, b) = (a, 0)$ satisfies the conditions. Otherwise, we have $\text{lc}(\tilde{g}) = a \star f_{d-1} + g_{d-2}$ since $a \star f_{d-1} \in L_{a \star f_d} = L_{gd-1}$ and $g_{d-2} \in L_{gd-1}$ we deduce that $\text{lc}(\tilde{g}) \in L_{gd-1}$. From Thm. 2.3.2, there exists a unique polynomial $b \in \mathbb{K}[x]$ with $\text{supp}(b) \subset S_{gd-1}$ such that $\text{lc}(\tilde{g}) = b \star g_{d-1}$. Hence, by construction $r = a \cdot f + (y-b) \cdot g$ has degree $< \text{deg}_y(g)$. For the uniqueness of $(a, b)$, consider $(a', b')$ another pair, which gives $\text{deg}_y((a-a') \cdot f + (b'-b) \cdot g) < \text{deg}_y(g)$. So, $(a-a') \star f_d = 0$ and $a = a'$ by Thm. 2.3.2. Finally, we must have $(b-b') \star g_{d-1} = 0$, so $b = b'$ again by Thm. 2.3.2. $\square$

For $0 \le j < d_y$, we consider the reverse truncated formal power
series $S_j = \sum_{k=j}^{D_y} u_{*,k}^{(j)} y^{D_y-k}$ representing the sequence $u^{(j)}$ at
precision $D_y$ with $D_y \ge 2d_y$. Also, we note $S_{-1} = \mathbf{u}_{*,0}^{(0)} y^{D_y+1}$.

LEMMA 3.2.2. For $0 \le j < d_y$, the pseudo-Euclidean division of $S_{j-1}$ by $S_j$ is well-defined.

PROOF. For $1 \le j < d_y$, from Thms. 3.1.6 and 3.1.7 and the con-
struction of $S_{j-1}$ and $S_j$ we deduce that the hypotheses of Thm. 3.2.1
are satisfied. For $j=0$, by construction of $S_{-1}$ the hypotheses of
Thm. 3.2.1 are also satisfied. Hence for $0 \le j < d_y$, the pseudo-
Euclidean division of $S_{j-1}$ by $S_j$ is well-defined. $\square$

The remainder of the pseudo-Euclidean division of $S_{j-1}$ by $S_j$ is
not exactly $S_{j+1}$ but has the same leading terms.

Definition 3.2.3 ([13, §11.1]). For a polynomial $p = \sum_{j=0}^{d} p_j y^j \in \mathbb{K}^N[y]$ of degree $d$ in $y$ and $k \le d$, we note $p \upharpoonright_k = \sum_{j=0}^{k} p_{d-j} y^{k-j}$ and $p \upharpoonright_k = y^{k-d} p$ when $k > d$.

LEMMA 3.2.4. Let $k \ge 1$. For $g \in \mathbb{K}[x, y]$ with $\text{deg}_y(g) = d \le k \le \text{deg}_y(S_0) = D_y$ and $w = g \star u$, we have $g \cdot S_0 \upharpoonright_k = f_{<d} + \sum_{j=0}^{k-d} w_{*,j} y^{k-j} + f_{>k}$ with $\text{deg}_y(f_{<d}) < d$ and $y^{k+1}$ divides $f_{>k}$.

PROOF. Let $p = g \cdot S_0 \upharpoonright_k \in \mathbb{K}^N[y]$. If we note $g_{*,l}$ the polynomial in $\mathbb{K}[x]$ associated to the monomial $y^l$ then from the arithmetic on $\mathbb{K}^N[y]$ defined in §3.2, we have $p_{k-j} = \sum_{\ell=0}^{d} g_{*,\ell} \star u_{*,j+\ell} = (g \star u)_{*,j}$ for $0 \le j \le k-d$. $\square$

LEMMA 3.2.5. For $0 \le j < d_y - 1$, if the pseudo-Euclidean division of $S_{j-1}$ by $S_j$ is $(c_j, d_j, \tilde{S}_{j+1})$ then we have $c_j = \bar{a}_j, d_j = \bar{b}_j$ with $(\bar{a}_j, \bar{b}_j)$ defined in Thm. 3.1.7 and $\tilde{S}_{j+1} \upharpoonright_{D_y-(j+1)-1} = S_{j+1} \upharpoonright_{D_y-(j+1)-1}$.

PROOF. For $j \ne 0$, the leading terms of $S_{j-1}$ are $u_{*,j-1}^{(j-1)} y^{D_y-j+1} + u_{*,j-2}^{(j-1)} y^{D_y-j}$ and similarly for $S_j$ we deduce from the proof of Thm. 3.2.1 that $c_j = \bar{a}_j$ and $d_j = \bar{b}_j$. For $j=0$, we have $S_{-1} = u_{*,0}^{(j-1)} y^{D_y+1}$ and $S_0$ has leading terms $u_{*,0}y^{D_y} + u_{*,1}y^{D_y-1}$ so we deduce that $c_0 = \bar{a}_0$ and $d_0 = \bar{b}_0$.

For $1 \le j < d_y - 1$, we have on the one hand the relation
$$ \mathbf{u}_{*,k}^{(j+1)} = \bar{a}_j \star \mathbf{u}_{*,k}^{(j-1)} + \mathbf{u}_{*,k-1}^{(j)} - \bar{b}_j \star \mathbf{u}_{*,k}^{(j)} $$
for $k > 0$ from Thm. 3.1.7.
On the other hand, we have $\tilde{S}_{j+1} = \bar{a}_j S_{j-1} + (y - \bar{b}_j) S_j$ which gives
$$ \tilde{S}_{j+1} = (\bar{a}_j \star \mathbf{u}_{*,j-1}^{(j-1)} + \mathbf{u}_{*,j}^{(j)}) y^{D_y-j+1} + \sum_{k=j}^{D_y-1} (\bar{a}_j \star \mathbf{u}_{*,k}^{(j-1)} + \mathbf{u}_{*,k-1}^{(j)} - \bar{b}_j \star \mathbf{u}_{*,k}^{(j)}) y^{D_y-k} + (\bar{a}_j \star \mathbf{u}_{*,D_y}^{(j-1)} - \bar{b}_j \star \mathbf{u}_{*,D_y}^{(j)}) $$
By definition of $\bar{a}_j$ and from Thm. 3.1.7, we deduce that $\text{deg}_y(\tilde{S}_{j+1}) = D_y - j - 1$ thus
$\tilde{S}_{j+1} \upharpoonright_{D_y-(j+1)-1} = S_{j+1} \upharpoonright_{D_y-(j+1)-1}$. For $j=0$, we apply the same arguments and obtain $\tilde{S}_1 \upharpoonright_{D_y-2} = S_1 \upharpoonright_{D_y-2}$. $\square$

For the purposes of Lms. 3.2.6 and 3.2.7, let $r_0, r_1, r'_0, r'_1 \in \mathbb{K}^N[y]$ and $k \ge 1$ such that $r_0 \upharpoonright_{2k} = r'_0 \upharpoonright_{2k}$ and $r_1 \upharpoonright_{2k-1} = r'_1 \upharpoonright_{2k-1}$. Assume that $d := \text{deg}_y(r_0) = \text{deg}_y(r_1)+1$ and $d' := \text{deg}_y(r'_0) = \text{deg}_y(r'_1)+1$.

LEMMA 3.2.6. Suppose that the pseudo-Euclidean division $(a_1, b_1, r_2)$ of $r_0$ by $r_1$ is well-defined, and that $\text{deg}_y(r_2) = d-2$. Then, the pseudo-Euclidean division $(a'_1, b'_1, r'_2)$ of $r'_0$ by $r'_1$ is also well-defined, and satisfies $a_1 = a'_1, b_1 = b'_1$, and $r_2 \upharpoonright_{2(k-1)-1} = r'_2 \upharpoonright_{2(k-1)-1}$. Moreover, $\text{deg}_y(r'_2) = d' - 2$ provided that $k \ge 2$.

PROOF. By assumption, the two leading terms of $r_0$ and $r'_0$ match,
and the same for $r_1$ and $r'_1$. Yet, the conditions (i), (ii), (iii) of
Thm. 3.2.1 which determine if a pseudo-Euclidean division is well-
defined only depends on the two leading terms of the dividend and
the divisor. In fact, $a_1, b_1$ only depend on those same two leading
terms. As a consequence, the pseudo-Euclidean division of $r'_0$ by $r'_1$
is well-defined, and $a_1 = a'_1, b_1 = b'_1$.

Assume w.l.o.g. $d' \le d$. The hypothesis $r_0 \upharpoonright_{2k} = r'_0 \upharpoonright_{2k}$ can be rewritten as $\text{deg}_y(r_0 - r'_0 y^{d-d'}) \le d - 2k - 1$. Likewise, $\text{deg}_y(r_1 - r'_1 y^{d-d'}) \le d - 2k - 1$. Considering that $r_2 = a_1 r_0 + (y - b_1)r_1$ and similarly for $r'_2$, we obtain that $\text{deg}_y(r_2 - r'_2 y^{d-d'}) \le d - 2k$. Whenever $k \ge 2$, $\text{deg}_y(r_2) = d - 2 > d - 2k \ge \text{deg}_y(r_2 - r'_2 y^{d-d'})$, which can only happen when $\text{deg}_y(r_2) = \text{deg}_y(r'_2) + d - d'$, i.e. $\text{deg}_y(r'_2) = d' - 2$, and $r_2 \upharpoonright_{2(k-1)-1} = r'_2 \upharpoonright_{2(k-1)-1}$. $\square$

LEMMA 3.2.7. Suppose that the first $k$ pseudo-Euclidean divisions $(a_j, b_j, r_{j+1})_{1 \le j \le k}$ starting from $r_0$ and $r_1$ are well-defined, and that $\text{deg}_y(r_j) = d-j$ for $1 \le j \le k$.

Then the first $k$ pseudo-Euclidean divisions $(a'_j, b'_j, r'_{j+1})_{1 \le j \le k}$ starting from $r'_0$ and $r'_1$ are also well-defined, and $a_j = a'_j, b_j = b'_j$, $\text{deg}_y(r'_j) = d' - j$ for $1 \le j \le k$. Moreover, $r_{j+1} \upharpoonright_{2(k-j)-1} = r'_{j+1} \upharpoonright_{2(k-j)-1}$ for $1 \le j < k$.

PROOF. Let us prove this statement by induction on $k$. The base case $k=1$ is a direct consequence of Lm. 3.2.6. For the induction step, suppose that $k \ge 2$ and that the lemma holds for $k-1$. Lm. 3.2.6 states that the first pseudo-Euclidean divisions $(a'_1, b'_1, r'_2)$ starting from $r'_0$ and $r'_1$ is well-defined, $\text{deg}_y(r'_2) = d' - 2$, $a'_1 = a_1, b'_1 = b_1$, and $r_2 \upharpoonright_{2(k-1)-1} = r'_2 \upharpoonright_{2(k-1)-1}$. It remains to apply our induction hypothesis to $k-1$ and $r_1, r_2, r'_1, r'_2$ to conclude. $\square$
If the *k* pseudo-Euclidean divisions ($a_j, b_j, r_{j+1}$)$_{0 \le j < k}$ starting from $r_{-1}$ and $r_0$ are well-defined then we have for $0 \le j < k$ the matrix relations $\begin{bmatrix} r_j \\ r_{j+1} \end{bmatrix} = Q_j \begin{bmatrix} r_{j-1} \\ r_j \end{bmatrix}$ where $Q_j := \begin{bmatrix} 0 & 1 \\ a_j & y-b_j \end{bmatrix}$. Thus, $\begin{bmatrix} r_j \\ r_{j+1} \end{bmatrix} = Q_j \cdots Q_0 \begin{bmatrix} r_{j-1} \\ r_0 \end{bmatrix}$ and by defining $\begin{bmatrix} s_j \\ s_{j+1} \end{bmatrix} \begin{bmatrix} t_j \\ t_{j+1} \end{bmatrix} = Q_j \cdots Q_0$, we have $s_{j+r_{-1}} + t_{j+r_0} = r_j$ for $0 \le j \le k$.

**THEOREM 3.2.8.** Let $k \ge 1$, $2k-1 \le D_y$ and $r_{-1}, r_0 \in \mathbb{K}^N[y]$ such that $r_{-1} = S_{-1} \upharpoonright_{2k}$ and $r_0 = S_0 \upharpoonright_{2k-1}$. Then, for all $0 \le j < \min(k, d_y)-1$, the pseudo-Euclidean division $(a_j, b_j, r_{j+1})$ of $r_{j-1}$ by $r_j$ is well-defined, $a_j = \bar{a}_j$ and $b_j = \bar{b}_j$ defined in Thm. 3.1.7, and also $r_j \upharpoonright_{2(k-j-1)} = S_j \upharpoonright_{2(k-j-1)}$ and $r_{j+1} \upharpoonright_{2(k-j-1)-1} = S_{j+1} \upharpoonright_{2(k-j-1)-1}$ with $\deg_y(r_{j+1}) = \deg_y(r_0) - (j+1)$. When $j = d_y-1$ and $2k-1 \ge 2d_y$, the pseudo-Euclidean division $(a_{d_y-1}, b_{d_y-1}, r_{d_y})$ on $r_{d_y-2}$ by $r_{d_y-1}$ is well-defined and $\deg_y(r_{d_y}) < \deg_y(r_0) - d_y$.

**PROOF.** We prove by induction for $0 \le j < \min(k, d_y) - 1$ that the pseudo-Euclidean division $(a_j, b_j, r_{j+1})$ of $r_{j-1}$ by $r_j$ is well-defined and $r_j \upharpoonright_{2(k-j-1)} = S_j \upharpoonright_{2(k-j-1)}$ and $r_{j+1} \upharpoonright_{2(k-j-1)-1} = S_{j+1} \upharpoonright_{2(k-j-1)-1}$. For $j > 0$, we suppose that the statement is true at step $j-1$ and we prove that it holds at step $j$. For every $j$, we have that $r_{j-1} \upharpoonright_{2(k-j)} = S_{j-1} \upharpoonright_{2(k-j)}$ and $r_j \upharpoonright_{2(k-j-1)} = S_j \upharpoonright_{2(k-j-1)}$ also from Lm. 3.2.5 the pseudo-Euclidean division $(\bar{a}_j, \bar{b}_j, \bar{S}_{j+1})$ of $S_{j-1}$ and $S_j$ is well-defined. Since $j < k-1$, we have $k-j \ge 2$ also by construction $\deg_y(S_j) = \deg_y(S_{j-1}) - 1$ so we can apply Lm. 3.2.6 and get that the pseudo-Euclidean $(a_j, b_j, r_{j+1})$ division of $r_{j-1}$ by $r_j$ is well-defined. On the one hand from Lm. 3.2.6, we have $r_{j+1} \upharpoonright_{2(k-j-1)-1} = \bar{S}_{j+1} \upharpoonright_{2(k-j-1)-1}$. On the other hand from Lm. 3.2.5, we have the equality $\bar{S}_{j+1} \upharpoonright_{D_y-(j+1)-1} = S_{j+1} \upharpoonright_{D_y-(j+1)-1}$. Since $2k-1 \le D_y$ and $j \ge 0$, we have $2k-1-2j-2 \le D_y-j-2$ so we conclude that $r_{j+1} \upharpoonright_{2(k-j-1)-1} = S_{j+1} \upharpoonright_{2(k-j-1)-1}$. Also from Lm. 3.2.6, we get $a_j = \bar{a}_j$ and $b_j = \bar{b}_j$ and $\deg_y(r_{j+1}) = \deg_y(r_0) - (j+1)$.

When $j = d_y - 1$ and $2k - 1 \ge 2d_y$ and the pseudo-Euclidean division $(\bar{a}_{d_y-1}, \bar{b}_{d_y-1}, r_{d_y})$ of $r_{d_y-2}$ by $r_{d_y-1}$ is well-defined by the same arguments so we have the relation $\bar{s}_{d_y}r_{j-1} + \bar{t}_{d_y}r_0 = r_{d_y}$. By hypothesis, we have $r_0 = S_0 \upharpoonright_{2k}$ so from Lm. 3.2.4 we can rewrite $\bar{t}_{d_y}r_0 = f_{<d_y} + \sum_{j=0}^{2k-1-d_y} \mathbf{u}_{*,j}^{(d_y)} y^{2k-1-j} + f_{>2k-1}$ with $\deg_y(f_{<d_y}) < d_y$ and $y^{2k}$ divides $f_{>2k-1}$. We deduce from the division property that $\deg_y(r_{d_y}) \le \deg_y(r_{d_y-1}) - 1 = 2k-1-d_y$ so by identification on the monomial basis we deduce that $\bar{s}_{j}r_{j-1} = f_{>2k-1}$ also since $\mathbf{u}^{(d_y)} = \mathbf{0}$ we have $\deg_y(r_{d_y}) < \deg_y(f_{<d_y}) < d_y$. Since $2k-1 \ge 2d_y$, it implies that $d_y \le 2k-1-d_y = \deg_y(r_0) - d_y$ hence $\deg_y(r_{d_y}) < \deg_y(r_0) - d_y$. $\square$

### 3.3 From successive remainders to C-relations

Let $r_{-1} = S_{-1}$ and $r_0 = S_0$. The definition of $S_{-1}$ is motivated by Thm. 3.2.1 and Lm. 3.2.4. Consider the successive remainders $(r_{-1}, r_0, \dots, r_{d_y})$ and relations $r_j = s_{j}r_{j-1} + t_j r_0$.

**LEMMA 3.3.1.** For $0 \le j < d_y$, we have $I(\mathrm{lc}(r_j)) = (I(\mathbf{u}) : \langle t_j \rangle) \cap \mathbb{K}[x]$ and $\langle 1 \rangle = \mathbb{K}[x] = (I(\mathbf{u}) : \langle t_{d_y} \rangle) \cap \mathbb{K}[x]$.

**PROOF.** Let $0 \le j < d_y$, from Lm. 3.1.8 and Thm. 3.2.8 we have $\mathbf{u}^{(j)} = t_j \star \mathbf{u}$ and $\mathrm{lc}(r_j) = \mathbf{u}_{*,j}^{(j)}$. From Thm. 3.1.6, we deduce that

$I(\mathrm{lc}(r_j)) = I(\mathbf{u}_{*,j}^{(j)}) = I(\mathbf{u}^{(j)}) \cap \mathbb{K}[x] = (I(\mathbf{u}) : \langle t_j \rangle) \cap \mathbb{K}[x]$. Also from Lm. 3.1.8, since $t_{d_y} \in I(\mathbf{u})$ we deduce the equality. $\square$

For $0 \le j < d_y$, we note $f_j$ be s.t. $\langle f_j \rangle = I(\mathrm{lc}(r_j))$ and $f_{d_y} = 1$.

**THEOREM 3.3.2.** The set $\{f_j t_j\}_{0 \le j \le d_y}$ is a Gröbner basis of $I(\mathbf{u})$.

**PROOF.** We verify that $f_j t_j \in I(\mathbf{u})$ and $\langle \mathrm{lm}(f_j t_j) \rangle_j = \langle \mathrm{lm}(I(\mathbf{u})) \rangle$. For $0 \le j \le d_y$, from Lm. 3.3.1 the polynomial $f_j t_j \in I(\mathbf{u})$. For $x^r y^s \in \mathrm{lm}(I(\mathbf{u}))$, if $s \ge d_y$ then $x^r y^s = \mathrm{lm}(x^r y^s - d_y t_{d_y})$. Otherwise, if $s < d_y$, by Thm. 3.1.5 we can find $f_s \in \mathbb{K}[x]$ such that $f_s p_{y^s} \in I(\mathbf{u})$ and $\mathrm{lm}(f_s p_{y^s}) = x^r y^s$. Since $p_{y^s} = t_s \mathrm{rem}(\mathcal{G}\mathbf{u})$ by Thm. 3.2.8, we deduce that $f_s t_s \in I(\mathbf{u})$. Note that since $\mathrm{lm}(t_s) = y^s = \mathrm{lm}(p_{y^s})$, $f_s t_s$ still has leading term $x^r y^s$. $\square$

From a Gröbner basis of $I(\mathbf{u})$, one can compute a minimal Gröbner basis of $I(\mathbf{u})$ with the following corollary.

**COROLLARY 3.3.3.** For $1 \le j \le d_y$, either $\mathrm{lm}(f_j t_j) \in \mathrm{lm}(\mathcal{G}\mathbf{u})$ or $\mathrm{deg}(f_j) = \mathrm{deg}(f_j)$.

**PROOF.** From the definition of minimal Gröbner basis, if $l \neq j$ and $\mathrm{lm}(f_l t_l)$ divides $\mathrm{lm}(f_j t_j)$ then $l \le j$ and $\mathrm{deg}(f_l) \le \mathrm{deg}(f_j)$. If $(\mathrm{deg}(f_j))_j$ is a decreasing sequence then it proves the claim.

First, we prove that $\mathrm{deg}(f_j) = \min\{\{r \mid x^r y^j \in \mathrm{lm}(I(\mathbf{u}))\}\}$. By definition $\langle f_j \rangle = I(\mathrm{lc}(r_j)) = (I(\mathbf{u}) : \langle t_j \rangle) \cap \mathbb{K}[x]$ by Lm. 3.3.1 so $\langle f_j \rangle = I(\mathbf{u}^{(j)}) \cap \mathbb{K}[x]$ from Thm. 2.3.2 and $\mathbf{u}^{(j)} = t_j \star \mathbf{u}$. Finally by Thm. 3.1.5, we can deduce that $I(\mathbf{u}^{(j)}) \cap \mathbb{K}[x] = \{r \mid x^r y^j \in \mathrm{lm}(I(\mathbf{u}))\}$. To conclude, if $x^{d_j} = \mathrm{deg}(f_j)$ then we have $\mathrm{lm}(f_j t_j) = (x^{d_j} y^j) y = x^{d_j} y^{j+1} \in \mathrm{lm}(I(\mathbf{u}))$ so $\mathrm{deg}(f_{j+1}) \le \mathrm{deg}(f_j)$. $\square$

## 4 ALGORITHMS

In the previous sections, we have considered bi-indexed sequences either as plain sequences $\mathbf{u} = (u_{i,j})_{i,j \in \mathbb{N}} \in \mathbb{K}^{N^2}$, or as polynomials with sequence coefficients $\mathbb{K}^N[y]$ in order to get relations out of a pseudo-Euclidean algorithm. At the moment, with the aim of fully describing our algorithms, we need to specify how the operations in $\mathbb{K}^N[y]$ are to be performed. Finite exact representations of uni-indexed sequences include the representation by the initial $d_x$ terms and the minimal relation, or the representation with the first $D_x \ge 2d_x$ terms, so that we can recover the relation. We choose the latter representation, and map these first $D_x$ terms in a reverse truncated formal power series as in [5, 6]. Doing so, the action $t \star \mathbf{u}$ can be computed using bivariate polynomial multiplication, which allows us to design efficient algorithms.

### 4.1 A finite polynomial representation

Let $\mathbf{u} = (u_{i,j})_{i,j \in \mathbb{N}}$ be a C-recursive sequence s.t. $x^{d_x}, y^{d_y} \in \mathrm{lm}(\mathcal{G}\mathbf{u})$ and consider bounds $D_x \ge 2d_x$ and $D_y \ge 2d_y$.

**Definition 4.1.1.** Fix $D_x \ge 2d_x$ and $D_y \ge 2d_y$. A polynomial $r \in \mathbb{K}[x, y]$ represents $\mathbf{u}$ at precision $(d, \delta)$ if $r \in \mathbb{K}[x, y] \le (D_x, D_y)$ and $r_{D_x-i, D_y-j} = u_{i,j}$ for $0 \le i \le d$ and $0 \le j \le \delta$.

For a representation $q$ of $\nu$ at precision $(D_x, \Delta)$, the addition term by term gives $q + r$, a representation of $\mathbf{u} + \nu$ at precision $(D_x, \min(\delta, \Delta))$. However, for the multiplication by a polynomial in $\mathbb{K}[x, y]$, we have to handle the same problem as in $\mathbb{K}^N[y]$ described in Lm. 3.2.4.

**LEMMA 4.1.2.** Let $r$ be a representation of $\mathbf{u}$ at precision $(d, \delta)$ and $t \in \mathbb{K}[x, y] \le (e, f)$. The polynomial $p = rt \mathrm{rem}\{x^{D_x+1}, y^{D_y+1}\}$ is a representation of $t \star \mathbf{u}$ at precision $(d - e, \delta - f)$.
**PROOF.** For $0 \le i \le d - e$ and $0 \le j \le \delta - f$, we have

$$p_{D_x-i, D_y-j} = \sum_{0 \le k \le e} \sum_{0 \le \ell \le f} t_{k, \ell} r_{D_x-i-k, D_y-j-\ell}$$

and $r_{D_x-i-k, D_y-j-\ell} = u_{i+k, j+\ell}$ by definition of a representation, so $p_{D_x-i, D_y-j} = (t \star u)_{i,j}$. By construction, $p \in \mathbb{K}[x, y] \le (D_x, D_y)$. $\square$

To handle the problem of decreasing precision in $x$, we can use fast univariate arithmetic to recover precision in $x$ when the C-relation $f_0 \in \mathcal{G}_u \cap \mathbb{K}[x]$ is known.

**THEOREM 4.1.3.** Let $r$ be a representation of $u$ at precision $(D_x, \delta)$ and $t \in \mathbb{K}[x, y] \le (e, f)$ with $e \le D_x$ and $f \le \delta$. From $r, t$ and $f_0$, one can compute a representation $p$, also denoted $t \cdot f_0 r$, of $w := t \star u$ at precision $(D_x, \delta - f)$ in $\tilde{O}(D_x \delta)$ operations in $\mathbb{K}$.

**PROOF.** Let $\tilde{t} = t \text{ rem } f_0$ with $\text{deg}_x(\tilde{t}) < d_x$. Since $f_0 \in I(u)$, we deduce that $w = \tilde{t} \star u$. The reduction of $t$ by $f_0$ requires $\tilde{O}(D_x \delta)$ operations by Lm. 2.1.1. Then, computing a representation $p$ of $w$ at precision $(D_x - d_x, \delta - f)$ has complexity in $\tilde{O}(D_x \delta)$ using $\tilde{t}$ and $r$. Since $D_x - d_x \ge d_x$, we can extend $p$ with $f_0$ from the univariate extension of Lm. 2.1.1 in $\tilde{O}(D_x \delta)$ operations. $\square$

We extend this definition of $f_0$ to matrix-vector multiplication with entries in $\mathbb{K}[x, y]$. With this new operation on representations in $\mathbb{K}[x, y]$, we can mimic the $\mathbb{K}[x, y]$ action on $\mathbb{K}^N[y]$ and apply a pseudo-Euclidean algorithm to solve the guessing problem on $u$.

## 4.2 Quotient algorithm

We now define the quotient algorithm of our pseudo-Euclidean division. For that, we need two subroutines for C-recursive unindexed sequences. The first one, GUESSINGUNIVAR($r$) takes a representation $r \in \mathbb{K}[x]$ of a C-recursive sequence $u$ at precision $D_x$ with $D_x \ge 2d_x$ and outputs the C-relation $f \in \mathbb{K}[x] \le d_x$ satisfying $\langle f \rangle = I(u)$. The other one, HANKELSOLVER($q, r, f$), takes a representation $q \in \mathbb{K}[x]$ of a C-recursive sequence $u$ at precision $D_x \ge 2(d_x - 1)$; a representation $r \in \mathbb{K}[x]$ of $v \in L_u$ at precision $d \ge d_x - 1$ and a C-relation $f \in \mathbb{K}[x] \le d_x$ s.t. $\langle f \rangle = I(u)$, and outputs $b \in \mathbb{K}[x] \le d_x$ satisfying $b \star u = v$. Both subroutines have complexities $\tilde{O}(D_x)$ (see §2.1).

**Algorithm 1 QUOBIVAR($f, g$)**

**Input:** Polynomials $f = \sum_{j=0}^{d} f_j(x)y^j$ and $g = \sum_{j=0}^{d-1} g_j(x)y^j$ satisfying the hypotheses of Thm. 3.2.1 when viewed in $\mathbb{K}^N[y]$ using the representation of $D_x + 1$ initial terms.

**Output:** $Q \in \mathbb{K}[x, y]^{2 \times 2}, \{p_1\} \subset \mathbb{K}[x]$ be s.t. $[\begin{matrix} g \\ r \end{matrix}] = Q \cdot p_0 \cdot [\begin{matrix} g \\ f \end{matrix}]$ with $\text{deg}_y(r) < \text{deg}_y(g)$ and $\langle p_1 \rangle = I(g_{d-1})$.

1. $p_0 \leftarrow \text{GUESSINGUNIVAR}(f_d(x))$

2. $a \leftarrow \text{HANKELSOLVER}(-f_d(x), g_{d-1}(x), f_0)$

3. $h(x) \leftarrow a \cdot p_0 \cdot f_{d-1}(x) + g_{d-2}(x)$

4. $p_1 \leftarrow \text{GUESSINGUNIVAR}(g_{d-1}(x))$

5. $b \leftarrow \text{HANKELSOLVER}(g_{d-1}(x), h(x), p_1)$

6. **return** $[\begin{matrix} 0 & 1 \\ a & y-b \end{matrix}]$, $\{p_1\}$

**LEMMA 4.2.1.** *QUOBIVAR is correct and has complexity in $\tilde{O}(D_x)$.*

**PROOF.** The polynomials $f, g$ viewed in $\mathbb{K}^N[y]$ satisfy the hypotheses of Thm. 3.2.1 so we consider $f_d, f_{d-1}, g_{d-1}, g_{d-2} \in \mathbb{K}^N$ the sequences represented by $f_d(x), f_{d-1}(x), g_{d}(x), g_{d-1}(x) \in \mathbb{K}[x]$.

Thm. 3.2.1 shows that there exists $a \in \mathbb{K}[x]_{<\deg(p_0)}$ such that $-a \star f_d = g_{d-1}$, that the call to HANKELSOLVER($-f_d(x), g_{d-1}(x), p_0$) computes. The update polynomial $h(x)$ represents the sequence $a \star f_{d-1} + g_{d-2}$. From Thm. 3.2.1, we can compute $b \in \mathbb{K}[x]_{<\deg(p_1)}$ such that $b \star g_{d-1} = a \star f_{d-1} + g_{d-2}$ also computed by the call to HANKELSOLVER($g_{d-1}(x), h(x), p_1$). By hypothesis of Thm. 3.2.1, $p_0$ is a C-relation on $f_d, f_{d-1}, g_{d-1}, g_{d-2}$, which ensures that $\text{deg}_y(r) < \text{deg}_y(g)$ by construction of the quotient matrix $Q$. Also, $\langle p_1 \rangle = I(g_{d-1})$ from the correctness of GUESSINGUNIVAR.

Computing $p_0, p_1 \in \mathbb{K}[x] \le d_x$ and $a, b \in \mathbb{K}[x] \le d_x$ have complexity in $\tilde{O}(D_x)$. The computation of $h(x)$ corresponds to univariate polynomial multiplication and addition of degree at most $D_x$ so it requires $\tilde{O}(D_x)$. Hence, we can bound the complexity of QUOBIVAR($f, g$) in $\tilde{O}(D_x)$. $\square$

## 4.3 Recursive pseudo-Euclidean algorithm

Based on the half-gcd algorithm, we build a divide and conquer pseudo-Euclidean algorithm, following the exposition of [13, Alg. 11.4]. Since our pseudo-Euclidean division has specific hypotheses, we define a requirement on the input of our algorithm.

*Input Requirement.* For the input $(r_{-1}, r_0, f_0, k)$, $f_0 \in \mathbb{K}[x]$ is a C-relation on the sequences represented by $r_{-1}$ and $r_0$, and there exists $0 \le \ell \le k$ such that the $\ell$ first pseudo-Euclidean division of $r_{-1}$ by $r_0$ are well-defined, and $\text{deg}_y(r_{\ell-1}) - 1 > \text{deg}_y(r_\ell)$ if $\ell < k$.

This requirement is satisfied at each recursive call for the input $r_{-1} = \sum_{i=0}^{D_x} u_{i,0} x^{D_x-i} y^{D_y+1}$, $r_0 = \sum_{j=0}^{D_y} \sum_{i=0}^{D_x} u_{i,j} x^{D_x-i} y^{D_y-j}$, $f_0 \in I(u) \cap \mathbb{K}[x]$ and $k \in \mathbb{N}$, but is not necessary satisfied by random bivariate polynomials.

**Algorithm 2 HALF-GCD-SEQ($r_{-1}, r_0, f_0, k$)**

**Input:** Representations $r_{-1}, r_0 \in \mathbb{K}[x, y]$, a C-relation $f_0 \in \mathbb{K}[x]$ and $k \in \mathbb{N}$ satisfying the input requirement §4.3.

**Output:** $R \in \mathbb{K}[x, y]^{2 \times 2}$ s.t. $[\begin{matrix} r_{\ell-1} \\ r_\ell \end{matrix}] = R \cdot f_0 \cdot [\begin{matrix} r_{\ell-1} \\ r_0 \end{matrix}]$, $T = [Q_0, \dots, Q_{\ell-1}] \in (\mathbb{K}[x, y]^{2 \times 2})^\ell$ s.t. $R = Q_{\ell-1} \cdots Q_0 \text{ rem } f_0$ and $\mathcal{F} = [f_0, \dots, f_{\ell-1}] \subset \mathbb{K}[x]$ s.t. $\langle f_j \rangle = I(\mathbf{w}_{*,j}^{(j)})$ with $\mathbf{w}^{(j)}$ corresponds to the sequence represented by $r_j$.

1. **if** $k = 0$ **then return** $[\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}]$, $\{\}$, $\{\}$

2. $d \leftarrow [k/2]$, $d^* \leftarrow k - d$

3. $R, T, \mathcal{F}_0 \leftarrow \text{HALF-GCD-SEQ}(r_{-1} \uparrow 2(d-1), r_0 \uparrow 2(d-1)-1, f_0, d - 1)$

4. $[\begin{matrix} r_{d-2} \\ r_{d-1} \end{matrix}] \leftarrow R \cdot f_0 \cdot [\begin{matrix} r_{d-1} \\ r_0 \end{matrix}]$

5. **if** $\text{deg}_y(r_{d-2}) - 1 > \text{deg}_y(r_{d-1})$ **then return** $R, T, \mathcal{F}_0$

6. $Q_{d-1}, \{f_{d-1}\} \leftarrow \text{QUOBIVAR}(r_{d-2} \uparrow 2, r_{d-1} \uparrow 1)$

7. $[\begin{matrix} r_{d-1} \\ r_d \end{matrix}] \leftarrow Q_{d-1} \cdot f_0 \cdot [\begin{matrix} r_{d-2} \\ r_{d-1} \end{matrix}]$

8. $S, U, \mathcal{F}_1 \leftarrow \text{HALF-GCD-SEQ}(r_{d-1} \uparrow 2d^*, r_d \uparrow 2d^*-1, f_0, d^*)$

9. **return** $(SQ_{d-1}R) \text{ rem } f_0, [T, Q_{d-1}, U], [\mathcal{F}_0, f_{d-1}, \mathcal{F}_1]$

We unroll the algorithm on a given sequence in [7].

**THEOREM 4.3.1.** *HALF-GCD-SEQ is correct. If $D_x$ (resp. $D_y + 1$) is the maximum degree in $x$ (resp. $y$) of $r_{-1}, r_0$ and $\lfloor D_y/2 \rfloor \le k \le D_y$ then HALF-GCD-SEQ($r_{-1}, r_0, f_0, k$) requires $\tilde{O}(D_x D_y)$ operations in $\mathbb{K}$.*

**PROOF.** We prove by induction on $j$, for any input $(r_{-1}, r_0, f_0, j)$ satisfying §4.3, HALF-GCD-SEQ is correct. For any input $(r_{-1}, r_0, f_0, 0)$ satisfying §4.3, the algorithm outputs $([\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}])$, $\{\}$, $\{\}$) which satisfies all the conditions of the algorithm output.
For $j \in \mathbb{N}$, we suppose the induction hypothesis at each step $i < j$ and we prove that the algorithm is correct for the input $(r_{-1}, r_0, f_0, j)$ satisfying §4.3. Consider the first $\ell$ pseudo-Euclidean divisions $(a_i, b_i, r_{i+1})_{0 \le i \le \ell-1}$ of $r_{-1}$ by $r_0$ with $0 \le \ell \le j$. Lm. 3.2.7 ensures that $(r_{-1} \upharpoonright_2 (d_{-1}), r_0 \upharpoonright_2 (d_{-1}-1), f_0, d-1)$ satisfies §4.3 so by the induction hypothesis we have the same quotient matrices $Q_i = \begin{bmatrix} 0 & 1 \\ a_i & y-b_i \end{bmatrix}$ for $0 \le i < \min(d-1, \ell)$. Since, $f_0$ is a C-relation on the sequences represented by $r_{-1}$ and $r_0$, it is a C-relation on $r_i$ due to $r_i = s_i r_{-1} + t_i r_0$ for $0 \le i < \min(d, \ell+1)$. If $\ell \le d-1$ then $R \cdot f_0 = \begin{bmatrix} r_{-1} \\ r_0 \end{bmatrix} = \begin{bmatrix} r_{\ell-1} \\ r_\ell \end{bmatrix}$ at Step 4 and $(R, T, \mathcal{F}_0)$ is the correct output. Otherwise, $r_{d-2}$ and $r_{d-1}$ are correctly computed at Step 4 from $r_{-1}, r_0$ and $R$ at precision $D_x$ in $x$. So, we can compute the quotient matrix $Q_d$ from $r_{d-2} \upharpoonright_2$ and $r_{d-1} \upharpoonright_1$ since the quotient algorithm only need the first two leading terms of each polynomial. The computation of $r_d$ from $r_{d-2}, r_{d-1}$ and $f_0$ is computed at full precision in $x$. From Lm. 3.2.7, since $j \ge \ell$, we have $d^* \ge \ell-d$, so the $(\ell-d)$ first pseudo-Euclidean division of $r_{d-1} \upharpoonright_2 d^*$ by $r_d \upharpoonright_2 d^*$-1 give the same results as the ones of $r_{d-1}$ and $r_d$. The second recursive call gives $S = Q_{\ell-1} \cdots Q_d$ rem $f_0, U = [Q_d, \dots, Q_{\ell-1}]$ and $\mathcal{F} = [f_d, \dots, f_{\ell-1}]$. Therefore, HALF-GCD-SEQ is correct.

For the complexity analysis, we suppose that $k$ is a power of 2 and we note $C(k)$ the cost of the computation. The base case requires $O(1)$ operations in $\mathbb{K}$. In the other cases due to the condition $k \in \Theta(D_y)$, the cost of the matrix multiplication $\cdot f_0$ is in $\tilde{O}(D_x k)$ and the call to QUOBIVAR is in $\tilde{O}(D_x)$ by Lm. 4.2.1. Finally, since the quotient matrices have all degree 1 in $y$, we deduce that $deg_y(S)$ and $deg_y(R)$ are less or equal to $\lfloor k/2 \rfloor$ and the degree in $x$ is bounded by $\tilde{O}(D_x)$ since the matrices are reduced by the relation $f_0$. Hence, the last matrix multiplication is in $\tilde{O}(D_x k)$. Note that the recursive calls continue to satisfy the condition $\lfloor D_y/2 \rfloor \le k \le D_y$. Thus, the cost $C(k)$ follows the recurrence $C(k) = 2C(k/2) + \tilde{O}(D_x k)$ and by the Master theorem we conclude that $C(k) = \tilde{O}(D_x D_y)$. $\square$

## 4.4 Guessing of bi-indexed sequences

From the result of Thm. 4.3.1, it remains to compute the cofactors $t_j$ from the quotient matrices to obtain a Gröbner basis of $I(\mathbf{u})$.

The product of matrices can be done recursively, we define RECURSIVEMATRIXPRODUCT$(T, k, f_0)$ with $T = [Q_0, \dots, Q_{\ell-1}] \in \mathbb{K}[x, y]^{2 \times 2}$ s.t. $deg_y(Q_j) = 1, 0 \le k < \ell$ and $f_0 \in \mathbb{K}[x]_{\le d_x}$ which computes the matrix $R = Q_k \cdots Q_0$ rem $f_0$. A call to RECURSIVEMATRI PRODUCT$(T, k, f_0)$ requires $\tilde{O}(d_x k)$ operations.

By combining the algorithms HALF-GCD-SEQ and RECURSIVE-MATRIXPRODUCT, we obtain a quasi-linear guessing algorithm for C-recursive bi-indexed sequences w.r.t. the lexicographic ordering.

The output of our algorithm is a minimal Gröbner basis with polynomials in $\mathbb{K}[x, y]_{\le (d_x, d_y)}$, so its size is in $O(|\mathcal{G}_u|d_x d_y)$.

**THEOREM 4.4.1.** *GUESSINGBIVAR is correct and has complexity in* $\tilde{O}(D_x D_y + |\mathcal{G}_u|d_x d_y)$.

**PROOF.** For the correctness, the polynomial $f_0$ computed from the call to GUESSINGONEVAR is in $\mathcal{G}_u \cap \mathbb{K}[x]$ since $I(u_{*,0}) = I(\mathbf{u}) \cap \mathbb{K}[x]$ by Thm. 3.1.6. From Thm. 3.2.8, we have that $(r_{-1}, r_0, f_0, k)$ satisfies §4.3 with $\ell = d_y$ and by Thm. 4.3.1, we deduce that $R = Q_{d_y-1} \cdots Q_0$ rem $f_0$, $T = [Q_0, \dots, Q_{d_y-1}]$ and $\mathcal{G}_x = [f_0, \dots, f_{d_y-1}]$ with $f_j \in I(\mathbf{u}_{*,j}^{(j)}) = I(\mathbf{u}^{(j)}) \cap \mathbb{K}[x]$ by Thm. 3.1.6. From Cor. 3.3.3, we only have to compute the relations which are not divisible by a

### Algorithm 3 GUESSINGBIVAR($\mathbf{u}$)

**Input:** The initial terms $(u_{i,j})_{0 \le i \le D_x, 0 \le j \le D_y}$ of a C-recursive sequence $\mathbf{u}$ satisfying Asm. A with $D_x \ge 2d_x$ and $D_y \ge 2d_y$.

**Output:** $\mathcal{G}$ a minimal Gröbner basis of $I(\mathbf{u})$ w.r.t. the order <.

1: $k \leftarrow [D_y/2]$

2: $r_{-1} \leftarrow \sum_{i=0}^{D_x} u_{i,0} x^{D_x-i} y^{D_y+1}, r_0 \leftarrow \sum_{j=0}^{D_y} \sum_{i=0}^{D_x} u_{i,j} x^{D_x-i} y^{D_y-j}$

3: $f_0 \leftarrow \text{GUESSINGUNIVAR}(\sum_{i=0}^{D_x} u_{i,0} x^{D_x-i})$

4: $R, T, \mathcal{G}_x \leftarrow \text{HALF-GCD-SEQ}(r_{-1}, r_0, f_0, k)$

5: $d \leftarrow \deg(f_0), \mathcal{G} \leftarrow \{f_0\}, j \leftarrow 1$

6: **for** $f \in \mathcal{G}_x \setminus \{f_0\}$ **do**

7:     **if** $\deg(f) < d$ **then**

8:         $\begin{bmatrix} s_j & t_j \\ s_{j+1} & t_{j+1} \end{bmatrix} \leftarrow \text{RECURSIVEMATRIXPRODUCT}(T, j, f_0)$

9:         $\mathcal{G} \leftarrow \mathcal{G} \cup \{(t_j f) \text{ rem } f_0\}, d \leftarrow \deg(f)$

10:     $j \leftarrow j + 1$

11:     $\mathcal{G} \leftarrow \mathcal{G} \cup \{t_{d_y}\}$

▷ $R = \begin{bmatrix} s_{d_y-1} & t_{d_y-1} \\ s_{d_y} & t_{d_y} \end{bmatrix}$ rem $f_0$

12:     **return** $\mathcal{G}$

previous one. For that, we distinguish them by the degree of $f_j$ and compute the corresponding cofactor $t_j$ when $\deg(f_j-1) \ne \deg(f_j)$. Finally, we get $t_{d_y}$ from the matrix $R$. Hence, GUESSINGBIVAR outputs a minimal Gröbner basis of $I(\mathbf{u})$ in $\mathcal{G}$.

For the complexity analysis, calling HALF-GCD-SEQ is in $\tilde{O}(D_x D_y)$ by Thm. 4.3.1. The loops on the polynomials of $\mathcal{G}_x$ add computations only if they compute a new polynomial in the minimal Gröbner basis and do at most $\tilde{O}(d_x d_y)$ operations. Finally, all the other instructions of the algorithm are in $\tilde{O}(D_x D_y)$. Hence, the complexity of GUESSINGBIVAR is in $\tilde{O}(D_x D_y + |\mathcal{G}_u|d_x d_y)$. $\square$

## 5 BENCHMARKS

The quasi-linearity of our guessing algorithm can be observed in practice from our implementation in MAPLE [7]. We compare the timings of our implementation also in MAPLE of guessing algorithms from [4, 19, 22]. For some we have to specialize the implementation for the lexicographic ordering with weighted degree ordering. For [4], we consider the adaptive version of the algorithms. We do not compare with [6], as under Asm. A, the computations are the same as in [4].

In our examples, we consider different shapes of staircase using Lazard's structure theorem [16] to build $\text{lm}(\mathcal{G}_u)$. We distinguish two particular shapes: simplex with $\text{lm}(\mathcal{G}_u) = \{x^{d_x-j} y^j\}_{0 \le j \le d_x}$ and L-shape with $\text{lm}(\mathcal{G}_u) = \{x^{d_x}, xy, y^{d_y}\}$. To begin with, we consider that we know $d_x, d_y$ and give exactly $D_x = 2d_x$ and $D_y = 2d_y$ in order to compute a minimal Gröbner basis of $I(\mathbf{u})$. The quantity size($\mathcal{G}_u$) corresponds to the number of coefficients in $\mathbb{K} = \mathbb{F}_{2^{16}+1}$ to represent $\mathcal{G}_u$. For simplex, Fig. 1 shows a quasi-linear growth on the timings of Alg. 3 following the growth of the quantity $|\mathcal{G}_u|d_x d_y$. For L-shape, the timings of Alg. 3 also follow the complexity found following the growth of the quantity $D_x D_y$. But it is outperformed by the adaptive version of the different algorithms. Next, for the second row of the Fig. 1 we now consider more initial terms of the sequence $\mathbf{u}$ than $D_x = 2d_x$ and $D_y = 2d_y$ by taking $(D_x, D_y) = (kd_x, kd_y)$ with $k \in \{10, 20, \dots, 90\}$. On Fig. 2 when $k > 80$, there is a crossover point on which the adaptive algorithm performs better, it is explained by the fact that these adaptive versions do not depend on the number of initial terms $D_x D_y$.
<table><caption>Figure 1: MAPLE implementation of several examples with initial terms ($D_x, D_y$) = ($2d_x, 2d_y$), timings in seconds.</caption><thead><tr><th></th><th><i>G<sub>u</sub></i>|<i>d<sub>x</sub></i>|<i>d<sub>y</sub></i></th><th><i>D<sub>x</sub></i>|<i>D<sub>y</sub></i></th><th>size(<i>G<sub>u</sub></i>)</th><th>[22]</th><th>[4]</th><th>[19]</th><th>Alg. 3</th></tr></thead><tbody><tr><td rowspan="7">simplex</td><td>27900</td><td>3600</td><td>9951</td><td>4452.3</td><td>20.5</td><td>204.4</td><td>4</td></tr><tr><td>127500</td><td>10000</td><td>44250</td><td>>24h</td><td>1394</td><td>3978.4</td><td>13.3</td></tr><tr><td>347900</td><td>19600</td><td>119348</td><td>>24h</td><td>17683.1</td><td>>10h</td><td>31</td></tr><tr><td>737100</td><td>32400</td><td>251248</td><td>>24h</td><td>>24h</td><td>>24h</td><td>62.2</td></tr><tr><td>1343100</td><td>48400</td><td>455937</td><td>>24h</td><td>>24h</td><td>>24h</td><td>115.5</td></tr><tr><td>2213900</td><td>67600</td><td>749439</td><td>>24h</td><td>>24h</td><td>>24h</td><td>175.5</td></tr><tr><td>3397500</td><td>90000</td><td>1147735</td><td>>24h</td><td>>24h</td><td>>24h</td><td>216.3</td></tr><tr><td rowspan="7">L-shape</td><td>4941900</td><td>115600</td><td>1666819</td><td>>24h</td><td>>24h</td><td>>24h</td><td>458.4</td></tr><tr><td>21600</td><td>28800</td><td>250</td><td>>24h</td><td>1</td><td>1425.8</td><td>14.4</td></tr><tr><td>117600</td><td>156800</td><td>570</td><td>>24h</td><td>27.5</td><td>>24h</td><td>234</td></tr><tr><td>290400</td><td>387200</td><td>890</td><td>>24h</td><td>114.7</td><td>>24h</td><td>465.5</td></tr><tr><td>540000</td><td>720000</td><td>1210</td><td>>24h</td><td>318.2</td><td>>24h</td><td>1021.3</td></tr><tr><td>866400</td><td>1155200</td><td>1530</td><td>>24h</td><td>749.1</td><td>>24h</td><td>1908</td></tr><tr><td>1269600</td><td>1692800</td><td>1850</td><td>>24h</td><td>1431.1</td><td>>24h</td><td>3583.1</td></tr></tbody></table>

<table><caption>Figure 2: MAPLE implementation of one example with initial terms ($D_x, D_y$) = ($kd_x, kd_y$), timings in seconds.</caption><thead><tr><th></th><th><i>G<sub>u</sub></i>|<i>d<sub>x</sub></i>|<i>d<sub>y</sub></i></th><th><i>D<sub>x</sub></i>|<i>D<sub>y</sub></i></th><th><i>k</i></th><th>[22]</th><th>[4]</th><th>[19]</th><th>Alg. 3</th></tr></thead><tbody><tr><td rowspan="8">simplex</td><td>127500</td><td>10000</td><td>2</td><td>>24h</td><td>1394</td><td>3978.4</td><td>13.3</td></tr><tr><td>127500</td><td>250000</td><td>10</td><td>>24h</td><td>1394</td><td>>24h</td><td>30.6</td></tr><tr><td>127500</td><td>1000000</td><td>20</td><td>>24h</td><td>1394</td><td>>24h</td><td>78</td></tr><tr><td>127500</td><td>2250000</td><td>30</td><td>>24h</td><td>1394</td><td>>24h</td><td>86.4</td></tr><tr><td>127500</td><td>4000000</td><td>40</td><td>>24h</td><td>1394</td><td>>24h</td><td>171.4</td></tr><tr><td>127500</td><td>6250000</td><td>50</td><td>>24h</td><td>1394</td><td>>24h</td><td>218.9</td></tr><tr><td>127500</td><td>9000000</td><td>60</td><td>>24h</td><td>1394</td><td>>24h</td><td>522.6</td></tr><tr><td>127500</td><td>12250000</td><td>70</td><td>>24h</td><td>1394</td><td>>24h</td><td>576.5</td></tr><tr><td>127500</td><td>16000000</td><td>80</td><td>>24h</td><td>1394</td><td>>24h</td><td>1008.1</td></tr><tr><td>127500</td><td>20250000</td><td>90</td><td>>24h</td><td>1394</td><td>>24h</td><td>2174.9</td></tr></tbody></table>

REFERENCES

[1] E. Becker, T. Mora, M. G. Marinari, and C. Traverso. 1994. The shape of the Shape Lemma. In Proceedings of the International Symposium on Symbolic and Algebraic Computation (Oxford, United Kingdom) (ISSAC '94). Association for Computing Machinery, New York, NY, USA, 129-133. https://doi.org/10.1145/190347.190382

[2] E. Berlekamp. 1968. Nonbinary BCH Decoding. IEEE Trans. Inform. Theory 14, 2 (1968), 242-242. https://doi.org/10.1109/TIT.1968.1054109

[3] J. Berthomieu, B. Boyer, and J.-Ch. Faugère. 2015. Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences. In Proceedings of the 2015 ACM on International Symposium on Symbolic and Algebraic Computation (Bath, United Kingdom) (ISSAC '15). ACM, New York, NY, USA, 61-68. https://doi.org/10.1145/2755996.2756673

[4] J. Berthomieu, B. Boyer, and J.-Ch. Faugère. 2017. Linear algebra for computing Gröbner bases of linear recursive multidimensional sequences. Journal of Symbolic Computation 83 (2017), 36-67. https://doi.org/10.1016/j.jsc.2016.11.005 Special issue on the conference ISSAC 2015.

[5] J. Berthomieu and J.-Ch. Faugère. 2018. A Polynomial-Division-Based Algorithm for Computing Linear Recurrence Relations. In Proceedings of the 2018 ACM International Symposium on Symbolic and Algebraic Computation (New York, NY, USA) (ISSAC '18). ACM, New York, NY, USA, 79-86. https://doi.org/10.1145/3208976.3209017

[6] J. Berthomieu and J.-Ch. Faugère. 2022. Polynomial-division-based algorithms for computing linear recurrence relations. Journal of Symbolic Computation 109 (2022), 1-30. https://doi.org/10.1016/j.jsc.2021.07.002

[7] J. Berthomieu, R. Lebreton, and K. Tran. 2025. CrecbiseqGuessing, a Maple Function for Quasi-Linear Guessing Minimal Lexicographic Gröbner bases of Ideals of C-Relations of Random Bi-Indexed Sequences. https://github.com/ktran11/CrecbiseqGuessing/blob/main/example/example_iterative.pdf.

[8] A. Bostan, G. Lecerf, and É. Schost. 2003. Tellegen's principle into practice. In Proceedings of the 2003 International Symposium on Symbolic and Algebraic Computation (Philadelphia, PA, USA) (ISSAC '03). Association for Computing Machinery, New York, NY, USA, 37-44. https://doi.org/10.1145/860854.860870

[9] R. P. Brent, F. G. Gustavson, and D. Y. Y. Yun. 1980. Fast solution of Toeplitz systems of equations and computation of Padé approximants. Journal of Algorithms 1, 3 (1980), 259-295. https://doi.org/10.1016/0196-6774(80)90013-9

[10] D. A. Cox, J. Little, and D. O'Shea. 2015. Ideals, Varieties, and Algorithms (4 ed.). Springer Cham. https://doi.org/10.1007/978-3-319-16721-3

[11] J-Ch. Faugère, P. Gianni, D. Lazard, and T. Mora. 1993. Efficient Computation of Zero-dimensional Gröbner Bases by Change of Ordering. Journal of Symbolic Computation 16, 4 (1993), 329-344. https://doi.org/10.1006/jSCO.1993.1051

[12] J.-Ch. Faugère and Ch. Mou. 2017. Sparse FGLM Algorithms. Journal of Symbolic Computation 80 (2017), 538-569. https://doi.org/10.1016/j.jsc.2016.07.025

[13] J. von zur Gathen and J. Gerhard. 2013. Modern Computer Algebra (3 ed.). Cambridge University Press. https://doi.org/10.1017/CBO9781139856065

[14] S. G. Hyun, V. Neiger, and É. Schost. 2021. Algorithms for Linearly Recurrent Sequences of Truncated Polynomials. In Proceedings of the 2021 International Symposium on Symbolic and Algebraic Computation (Virtual Event, Russian Federation) (ISSAC '21). Association for Computing Machinery, New York, NY, USA, 201-208. https://doi.org/10.1145/3452143.3465533

[15] D. E. Knuth. 1970. The analysis of algorithms. In Actes du Congres International des Mathématiciens (Nice, 1970), Vol. 3. 269-274.

[16] D. Lazard. 1985. Ideal Bases and Primary Decomposition: Case of Two Variables. J. Symb. Comput. 1, 3 (1985), 261-270. https://doi.org/10.1016/S0747-7171(85)80035-3

[17] J. Massey. 1969. Shift-register synthesis and BCH decoding. IEEE Trans. Inf. Theor. 15, 1 (1969), 122-127. https://doi.org/10.1109/TIT.1969.1054260

[18] R. T. Moenck. 1973. Fast Computation of GCDs. In Proceedings of the Fifth Annual ACM Symposium on Theory of Computing (Austin, Texas, USA) (STOC '73). Association for Computing Machinery, New York, NY, USA, 142-151. https://doi.org/10.1145/800125.804045

[19] B. Mourrain. 2017. Fast Algorithm for Border Bases of Artinian Gorenstein Algebras. In Proceedings of the 2017 ACM International Symposium on Symbolic and Algebraic Computation (Kaiserslautern, Germany) (ISSAC '17). Association for Computing Machinery, New York, NY, USA, 333-340. https://doi.org/10.1145/3087604.3087632

[20] S. Naldi and V. Neiger. 2020. A divide-and-conquer algorithm for computing Gröbner bases of syzygies in finite dimension. In Proceedings of the 45th International Symposium on Symbolic and Algebraic Computation (Kalamata, Greece) (ISSAC '20). Association for Computing Machinery, New York, NY, USA, 380-387. https://doi.org/10.1145/3373207.3404059

[21] Sh. Sakata. 1988. Finding a minimal set of linear recurring relations capable of generating a given finite two-dimensional array. J. Symbolic Comput. 5, 3 (1988), 321-337. https://doi.org/10.1016/S0747-7171(88)80033-6

[22] Sh. Sakata. 1990. Extension of the Berlekamp-Massey Algorithm to N Dimensions. Inform. and Comput. 84, 2 (1990), 207-239. https://doi.org/10.1016/0890-5401(90)90039-K

[23] Sh. Sakata. 2009. The BMS Algorithm. In Gröbner Bases, Coding, and Cryptography, Massimiliano Sala, Shojiro Sakata, Teo Mora, Carlo Traverso, and Ludovic Perret (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 143-163. https://doi.org/10.1007/978-3-540-93806-4_9

[24] A. Schönhage. 1971. Schnelle Berechnung von Kettenbruchentwicklungen. Acta Informatica 1, 2 (1971), 139-144. https://doi.org/10.1007/BF00289520

[25] G. Villard. 2018. On Computing the Resultant of Generic Bivariate Polynomials. In Proceedings of the 2018 ACM International Symposium on Symbolic and Algebraic Computation (New York, NY, USA) (ISSAC '18). Association for Computing Machinery, New York, NY, USA, 391-398. https://doi.org/10.1145/3208976.3209020

[26] G. Villard. 2025. Bivariate polynomial reduction and elimination ideal over finite fields. Journal of Symbolic Computation 127 (2025), 102367. https://doi.org/10.1016/j.jsc.2024.102367

[27] D. Wiedemann. 1986. Solving sparse linear equations over finite fields. IEEE Trans. Inf. Theory 32, 1 (1986), 54-62.