# Simple generators of rational function fields

Alexander Demin, Gleb Pogudin

► To cite this version:

Alexander Demin, Gleb Pogudin. Simple generators of rational function fields. 2026. hal-05505433v2

HAL Id: hal-05505433

https://hal.science/hal-05505433v2

Preprint submitted on 5 Mar 2026

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 CC BY 4.0 - Attribution - International License
# Simple generators of rational function fields

Alexander Demin,* Gleb Pogudin†

March 5, 2026

## Abstract

Consider a subfield of the field of rational functions in several indeterminates. We present an algorithm that, given a set of generators of such a subfield, finds a simple generating set. We provide an implementation of the algorithm and show that it improves upon the state of the art both in efficiency and the quality of the results. Furthermore, we demonstrate the utility of simplified generators through several case studies from different application domains, such as structural parameter identifiability. The main algorithmic novelties include performing only partial Gröbner basis computation via sparse interpolation and efficient search for polynomials of a fixed degree in a subfield of the rational function field.

## 1 Introduction

Let $k$ be a field, and consider the field of rational functions $k(\mathbf{x})$, where $\mathbf{x} = (x_1, \dots, x_n)$ are indeterminates. Intermediate subfields $k \subset \mathcal{E} \subset k(\mathbf{x})$ are ubiquitous in algebra and algebraic geometry. Their theoretical properties have been extensively studied during the 20-th century, in particular, due to the central role they play in the Lüroth problem [9] and Hilbert 14-th problem [83, 102].

A starting point for the computational study of the subfields of the rational function field is the fact that every such field can be generated over $k$ by finitely many elements [84, Ch. VIII, Ex. 4]; this contrasts sharply with the situation for intermediate subalgebras of the polynomial algebra $k \subset \mathcal{A} \subset k[\mathbf{x}]$ [112]. Such a finite generating set is then a natural data structure for representing subfields of $k(\mathbf{x})$. Within this framework, a number of algorithmic results have been obtained in the 90-s [77, 100, 101, 104, 121], in particular, a solution to the membership problem. One of the key technical tools introduced first by Ollivier [104, Ch. III, §2] and then by Müller-Quade and Steinwandt [100, 101] has been certain polynomial ideals in $k(\mathbf{x})[\mathbf{y}]$ with $\mathbf{y} = (y_1, \dots, y_n)$ which we call *OMS ideals* (Section 2.7); they play a prominent role in this paper as well. Another line of research has been devoted to subfields of transcendence degree one which have a particularly nice structure due to the Lüroth theorem [14, 54].

More recent advances in computing with rational function fields are often motivated by particular applications such as computer vision [13] or structural parameter identifiability of dynamical models [37, 105]. In the latter case, the field $\mathcal{E}$ of interest consists of quantities that may be inferred from the input-output data of the model. The existing algorithms for computing the generators of $\mathcal{E}$ typically produce extremely complicated expressions, which cannot be further analyzed and interpreted by a modeler. Therefore, one of the key questions in this context is the question of *generator simplification*, which can be stated as follows:

**Given:** $\mathbf{g} := (g_1, \dots, g_m) \in k(\mathbf{x})$, such that $\mathcal{E} = k(\mathbf{g})$;

**Find:** $\mathbf{h} := (h_1, \dots, h_\ell) \in k(\mathbf{x})$, such that $\mathcal{E} = k(\mathbf{h})$ and $\mathbf{h}$ are simpler than $\mathbf{g}$.

We will deliberately not give a formal definition of what being simple means, as we do not expect that there is a single mathematical definition capable of capturing the intuitive notions of simplicity used in applications. In general, it may be preferable to have a generating set of small cardinality, where each element is a sparse rational function of low degree with small coefficients. To give some intuition on how the simplification may look like and what benefits it may bring, let us give an illustrative example.

**Example 1.1.** Consider a triangle with side lengths $a, b, c$. Using Heron's formula, we can write the length of the altitude from $a$ as

$$h_a^2 = \frac{4s(s-a)(s-b)(s-c)}{a^2}, \quad \text{where } s = \frac{a+b+c}{2},$$

and analogous formulas hold for $h_b$ and $h_c$. By computing a simplified generating set of $\mathcal{E} := \mathbb{Q}(h_a^2, h_b^2, h_c^2)$, we find that $\mathcal{E} = \mathbb{Q}(a^2, b^2, c^2)$. Since the side lengths are always positive, we conclude that they can be determined from the altitude lengths confirming a well-known fact from elementary geometry.

<footnote>*National Research University Higher School of Economics, Moscow, Russia, asdemin_2@edu.hse.ru</footnote>

<footnote>†LIX, CNRS, École polytechnique, Institute Polytechnique de Paris, Paris, France, gleb.pogudin@polytechnique.edu</footnote>
Already the Lüroth problem can be viewed as a version of the simplification question focusing solely on the cardinality of the generating set. Specifically, any generating set of $k \subset \mathcal{E} \subset k(\mathbf{x})$ over $k$ must contain at least $\text{trdeg}_k \mathcal{E}$ elements. On the other hand, if $\text{char } k = 0$, the primitive element theorem implies that $\text{trdeg}_k \mathcal{E} + 1$ generators are always sufficient. Whether or not using only $\text{trdeg}_k \mathcal{E}$ generators is always possible is the Lüroth problem, and the answer is negative starting from $n = 3$. Furthermore, deciding whether $\text{trdeg}_k \mathcal{E} + 1$ generators are necessary is equivalent to determining if a unirational variety is rational, which is an old open problem in algebraic geometry [120]. For a related Noether problem, an approach based on searching for generators of small degree was used to give a computational solution in several interesting cases [77].

The question of simplification in its broader, more intuitive, form was already considered in [101] where it was proposed to use the coefficients of a reduced Gröbner basis of the OMS ideal of $\mathcal{E}$ as a “canonical generating set”. In the context of structural identifiability, this approach was further refined [105] by taking this canonical generating set and filtering out redundant generators. The resulting algorithm has been implemented in Maple [105], and the implementation was later incorporated into a web-application [67]. A construction reminiscent of OMS ideals was used in [92, 93, 94] to search for a simple transcendence basis of the relative algebraic closure of $\mathcal{E}$ inside $k(\mathbf{x})$ also in connection to the structural identifiability problem.

The main contribution of the present paper is a new algorithm that takes as input a set of generators of a subfield $\mathcal{E} \subset k(\mathbf{x})$ and returns a simplified set of generators, assuming $\text{char } k = 0$. We provide an open-source implementation of the algorithm for the case $k = \mathbb{Q}$ in Julia. Our approach also relies on OMS ideals but improves upon the state of the art both in terms of *efficiency* and *quality of simplification*. Let us describe our contributions along these two axes separately.

From the *efficiency* perspective, the key question is how to perform Gröbner basis computations with OMS ideals. These ideals belong to the ring $k(\mathbf{x})[\mathbf{y}]$ over the rational function field which makes computations with them expensive due to the intermediate expression swell. Our algorithm follows the evaluation-interpolation approach: it computes multiple Gröbner bases by specializing $\mathbf{x}$ and then reconstructs the coefficients by using sparse rational function interpolation [29, 124]. There are two rationale behind this choice of methodology.

First, interpolation is one of the standard ways to alleviate intermediate expression swell. It has been successfully employed for similar problems over the rational function field such as GCD computation [58, 124] and linear system solving [70] and explored for Gröbner basis computation [16, Ch. 5] and, recently, for generic parametrized polynomial system solving [27]. Our approach includes further optimizations such as computing over finite fields and using tracing to speed up repeated Gröbner basis computations (see Section 5.5). The resulting implementation is competitive for computing Gröbner bases in $\mathbb{Q}(\mathbf{x})[\mathbf{y}]$, significantly improving upon the prior interpolation-based algorithms, but may be slower than other, not interpolation-based methods (see Section 5.6).

And this is where the second key advantage of our method comes into play: it does not require computing the full Gröbner basis. Using evaluation-interpolation allows us to stop early, once only low-degree coefficients of the Gröbner basis have been interpolated, and these coefficients are the ones of interest in the simplification problem. So, we gradually increase the degree of interpolation stopping once the computed coefficients generate the whole field. Thus, compared to the earlier approach from [105], instead of computing all the coefficients and then filtering out the large redundant ones, we now avoid computing them at the first place. This strategy allows to reduce the number of evaluations and, consequently, the runtime dramatically and solve the problems which are out of reach otherwise (see Section 5.4).

We improve the *quality of simplification* by creating a pool of potential generators (which contains the coefficients of the OMS ideal) and employing a refined heuristic for selecting the final result from this pool. The main tool for adding extra candidates to the pool is our algorithm for finding polynomial elements of the subfield up to a fixed degree, which is of independent interest. While a theoretical algorithm for this task has appeared already in [101], it required computing the full Gröbner basis of the OMS ideal. In the spirit of the evaluation-based approach we follow, our algorithm performs only Gröbner basis computation on the specialized ideal which makes it efficient for practical computation.

Since we do not fix a formal notion of simplicity, assessing the *quality of simplification* is a delicate task. First, we run our algorithm and the one from [105] on a large suite of benchmarks (53 problems) coming from applications to structural identifiability and perform a manual inspection of results. We show that our algorithm is more efficient and argue that it produces simpler results (see Section 6.2). Furthermore, we perform a more detailed analysis of the output of our algorithm on several case studies coming from structural identifiability of ODE models, observability of discrete dynamics, and fields of invariants. In addition to syntactic simplicity, we discuss the usefulness of the simplification performed by our algorithm.

The rest of the paper is organized as follows. We review and adapt the main tools (e.g., sparse interpolation, OMS ideals) in Section 2. The main building blocks of our algorithm are presented in Section 3, and we use them to describe and justify our main algorithm in Section 4. Section 5 is devoted to implementation details including discussion of the impact of the main design choices on the performance. We analyze the
quality of simplification by our algorithm on a large benchmark suite in Section 6 (the details of the suite are summarized in Appendix A). We finish with a more in-depth discussion of a set of case studies coming from three different application domains in Section 7.

**Acknowledgements.** The authors would like to thank Hoon Hong, Joe Kileel, Grégoire Lecerf, Nicolette Meshkat, and François Ollivier for helpful discussions. The authors are also greatful to Jan Flusser and Evelyne Hubert for stimulating discussions of Section 7.3. GP was supported by the French ANR-22-CE48-0008 OCCAM and ANR-22-CE48-0016 NODE projects.

# 2 Preliminaries

## 2.1 Notation

Let $k$ denote an effective field, meaning that there exist algorithms for the arithmetic in $k$. For a field $k$, let $\bar{k}$ denote its algebraic closure. We will use shorthands for tuples, such as $\mathbf{x} := (x_1, \dots, x_n)$. For a multivariate polynomial $f \in k[\mathbf{x}]$, we will write $d_f := \deg f$ for the total degree of $f$ and $s_f$ for the number of nonzero terms in $f$. We define the degree of zero polynomial to be $-\infty$. We have $\mathbb{N}_0 := 0, 1, 2, \dots$.

Products, sums, and powering of tuples are performed component-wise, more precisely, for every $\gamma, \omega, \sigma \in k^n$ and every $\alpha \in \mathbb{N}_0^n$, we write $\gamma\omega^\alpha + \sigma := (\gamma_1\omega_1^{\alpha_1} + \sigma_1, \dots, \gamma_n\omega_n^{\alpha_n} + \sigma_n)$. We extend this notation to scalars by treating a scalar as a tuple of $n$ identical elements. In place of $\gamma, \omega, \sigma$ we may use indeterminates or tuples thereof. Given a tuple $\mathbf{a}$ (say, $\mathbf{a} \in k^n$), we write $\Pi(\mathbf{a}) := \prod_{i=1,\dots,n} a_i$.

The degrevlex monomial ordering is a total ordering on the set of monomials, which we may denote by $<_{\text{dr1}}$. For $\mathbf{x}^\alpha, \mathbf{x}^\beta$ with $\alpha, \beta \in \mathbb{N}_0^n$, we have $\mathbf{x}^\alpha <_{\text{dr1}} \mathbf{x}^\beta$ whenever the sum of elements of $\alpha$ is smaller than that of $\beta$, or, them being equal, the right-most nonzero entry of $\alpha - \beta$ is positive.

For polynomials $f_1, \dots, f_\ell$ from $k[\mathbf{x}]$, let $\langle f_1, \dots, f_\ell \rangle$ denote the ideal in $k[\mathbf{x}]$ generated by $f_1, \dots, f_\ell$. For an ideal $I$ of $k[\mathbf{x}]$ and an element $a \in k[\mathbf{x}]$, denote $I: a^\infty := \{r \in k[\mathbf{x}] \mid \exists \ell : ra^\ell \in I\}$. This set is also an ideal.

## 2.2 DeMillo-Lipton-Schwartz-Zippel lemma

Let $k$ be an effective field and $\mathbf{x} := (x_1, \dots, x_n)$ be an $n$-tuple of variables. Let the coordinates of $\gamma := (\gamma_1, \dots, \gamma_n)$ be chosen independently and uniformly from a finite subset $S \subseteq k$. Many of our probabilistic statements will rely on the following lemma:

**Lemma 2.1 (DeMillo-Lipton-Schwarz-Zippel [33, 117, 133]).** Let $f \in k[\mathbf{x}]$ be a nonzero polynomial of degree $d$. The probability that $f(\gamma) = 0$ is at most $d/|S|$.

**Corollary 2.2.** Let $f \in k[\mathbf{x}]$ be a nonzero polynomial of degree $d$. Let $\sigma_1, \dots, \sigma_T \in k^n$ be arbitrary but fixed. Let $\omega_1, \dots, \omega_T \in (k \setminus \{0\})^n$ be arbitrary but fixed. The probability that there exists $i = 1, \dots, T$ such that $f(\omega_i; \gamma + \sigma_i) = 0$ is at most $dT/|S|$.

*Proof.* Let $p := f(\omega_1\mathbf{x} + \sigma_1) \cdots f(\omega_T\mathbf{x} + \sigma_T)$. Note $p \neq 0$, since for every $i = 1, \dots, T$ none of the components of $\omega_i$ are zero. We apply Lemma 2.1 to $p$. $\square$

## 2.3 Univariate rational function interpolation

Let $F = A/B$ be a rational function in $k(x)$ with $A, B \in k[x]$, $\gcd(A, B) = 1$, and $d_A := \deg A$ and $d_B := \deg B$. Let $D_A \ge d_A$ and $D_B \ge d_B$ be the bounds on the total degrees of $A$ and $B$, respectively. Let $D$ be such that $D \ge D_A + D_B + 2$, and let $u_1, \dots, u_D \in k$ be pairwise different such that for every $i = 1, \dots, D$ we have $B(u_i) \neq 0$ (assuming $k$ is sufficiently large). Univariate rational function interpolation is the problem of recovering $F$ from the evaluations $F(u_1), \dots, F(u_D)$ and the bounds $D_A, D_B$.

The solution is unique up to a scaling of $A$ and $B$, and can be computed using univariate polynomial interpolation and the Extended Euclidean Algorithm [126, Corollary 5.18]. This method is also referred to as Cauchy interpolation.

**Remark 2.3.** When the bounds on the degrees $D_A$ and $D_B$ are not known, it is still possible to recover a solution using [79, Algorithm MQRFR] instead of the EEA. Later in this paper, these two methods may be used interchangeably.
## 2.4 Sparse polynomial interpolation

Let $f$ be a polynomial in $k[x]$ with $x := (x_1, \dots, x_n)$. Let $d_f$ be an integer, the total degree of $f$. Sparse interpolation is the problem of recovering the sparse representation of $f$ in the form

$$f = \sum_{i=1, \dots, s_f} a_i \Pi(x^{\alpha_i}) \qquad (1)$$

from a sufficient number of evaluations of $f$, where $s_f$ is the number of terms in $f$ and, for every $i = 1, \dots, s_f$, we have $a_i \in k \setminus \{0\}$ and $\alpha_i \in \mathbb{N}_0^n$ are pairwise different.

In the case an upper bound on the number of terms $s_f$ is known, the algorithm by Ben-Or and Tiwari [10] (which generalizes the Prony's method [108]) recovers a sparse representation of $f$ in the form (1) deterministically from the evaluations $f(\omega^0), f(\omega^1), \dots, f(\omega^{2T-1})$ with $T \ge s_f$ for some carefully chosen $\omega \in (k \setminus \{0\})^n$. We shall call such $\omega$ an *admissible ratio*, following [35]. A variant of the algorithm by Ben-Or and Tiwari is summarized in Algorithm 1.

There exist several strategies for choosing an admissible ratio $\omega$. A common requirement is that for every $\alpha \in \mathbb{N}_0^n$ such that the degree of $\Pi(x^\alpha)$ does not exceed $d_f$ it is possible to determine the exponent $\alpha$ from the product $\Pi(\omega^\alpha)$. In practice, when $k = \mathbb{Q}$ or when $k$ is a finite field of large characteristic, one may take $\omega := (\omega_1, \dots, \omega_n) = (2, 3, 5, \dots, p_n)$, where $p_n$ denotes the $n$-th prime number. In this case, the exponent can be recovered via integer factorization (more precisely, via trial division). Other approaches include reduction to univariate interpolation by Kronecker substitution [126, Section 8.4] and picking the evaluation point of a large order in a suitable extension of $k$ [125].

**Algorithm 1** Deterministic sparse polynomial interpolation [10, Figure 1], [125, Algorithm 4]

**Input:** an admissible ratio $\omega \in (k \setminus \{0\})^n$; the evaluations $f(\omega^0), f(\omega^1), \dots, f(\omega^{2T-1})$.

**Output:** if $T \ge s_f$, the sparse representation of $f$ as in (1); if $T < s_f$, an incorrect result.

(Step 1) Compute (via Páde approximation [126, Section 5.9]) the minimal $t \le T$, a monic $\Lambda \in k[z]$ of degree $t$, and $N \in k[z]$ of degree less than $t$, such that the following identity holds modulo $O(z^{2T})$:

$$\sum_{i=0}^{\infty} f(\omega^i) z^i = \frac{N(z)}{\Lambda(z)}.$$

(Step 2) Find the roots $\Pi(\omega^{-\alpha_1}), \dots, \Pi(\omega^{-\alpha_t})$ of $\Lambda(z)$ using polynomial root-finding.

(Step 3) Compute the exponents $\alpha_1, \dots, \alpha_t$ using integer factoring.

(Step 4) Compute the coefficients $a_1, \dots, a_t$ by solving the linear system (for example, using [74]):

$$a_1 + \dots + a_t = f(\omega^0)$$

...

$$a_1 \Pi(\omega^{(t-1)\alpha_1}) + \dots + a_t \Pi(\omega^{(t-1)\alpha_t}) = f(\omega^{t-1})$$

(Step 5) Return $\alpha_1, \dots, \alpha_t$ and $a_1, \dots, a_t$.

**Remark 2.4.** The fact that a field is effective does not immediately imply that there is an algorithm for finding root of univariate polynomials over it. Over the rationals and over finite fields there exist efficient deterministic algorithms [10, Figure 2], [21], [126, Section 14.5.], [53].

**Remark 2.5.** Computing an upper bound on the number of terms $s_f$ deterministically is a hard problem; see the discussion in [10, Section 8]. For this reason, in Algorithm 1, when $T < s_f$, we allow an incorrect result instead of a proper FAIL. An upper bound can be computed probabilistically [125, Section 4.6].

**Remark 2.6 (Diversification [48].)** In the following sections, we will be interested in estimating the probability that for a fixed nonzero $\tilde{f}(\mathbf{x}) := f(\mathbf{x} + \boldsymbol{\sigma})$ any of the evaluations $\tilde{f}(\omega_1), \dots, \tilde{f}(\omega_T)$ vanish at a fixed $\omega_1, \dots, \omega_T \in (k \setminus \{0\})^n$ and $\boldsymbol{\sigma} \in k^n$. In theoretical arguments, we may replace $\tilde{f}(\mathbf{x})$ with $\tilde{f}(\gamma\mathbf{x})$ and apply Corollary 2.2, where the coordinates of $\gamma \in k^n$ are chosen independently and uniformly from $S \subseteq k$. In practice, this amounts to evaluating $\tilde{f}(\gamma\omega_1), \dots, \tilde{f}(\gamma\omega_T)$; when the end goal is to interpolate the original $\tilde{f}(\mathbf{x})$ from these evaluations, we would additionally require that the coordinates of $\gamma$ are nonzero, so that the original polynomial could be restored by applying the change of variables $\mathbf{x} \mapsto \mathbf{x}/\gamma$.
## 2.5 Sparse rational function interpolation

Let $F$ be a rational function in $k(x)$ with $x := (x_1, \dots, x_n)$. Sparse interpolation is the problem of recovering the sparse representation of $F$ in the form

$$F = \frac{A}{B}, \quad A = \sum_{i=1, \dots, s_A} a_i \Pi(\mathbf{x}^{\alpha_i}), \quad B = \sum_{i=1, \dots, s_B} b_i \Pi(\mathbf{x}^{\beta_i}), \qquad (2)$$

from a sufficient number of evaluations of $F$. Here for every $i = 1, \dots, s_A$ we have $a_i \in k \setminus \{0\}$ and $\alpha_i \in \mathbb{N}_0^n$ are pairwise different and for every $j = 1, \dots, s_B$ we have $b_j \in k \setminus \{0\}$ and $\beta_j \in \mathbb{N}_0^n$ are pairwise different. We also assume that $\gcd(A, B) = 1$. Let $d_A$ and $d_B$ be the total degrees of $A$ and $B$, respectively.

Assume we have access to a blackbox that evaluates $F$ at a point in $k^n$. The algorithm by van der Hoeven and Lecerf [124, Section 3.3.] is a probabilistic algorithm for recovering the sparse representation of a rational function as in (2) from a sequence of evaluations of the blackbox. We will show that the probability of failure in this algorithm can be made arbitrarily close to zero assuming we can sample elements independently and uniformly from an arbitrarily large subset $S \subseteq k \setminus \{0\}$. For completeness, we describe the algorithm here.

Introduce a new indeterminate $x_0$. Consider the rational function

$$\hat{F}(x_0, \dots, x_n) := x_0^{d_A - d_B} \frac{A(x_1/x_0, \dots, x_n/x_0)}{B(x_1/x_0, \dots, x_n/x_0)},$$

where we also introduce for convenience

$$\hat{A} := x_0^{d_A} A(x_1/x_0, \dots, x_n/x_0),$$

$$\hat{B} := x_0^{d_B} B(x_1/x_0, \dots, x_n/x_0).$$

Note that $\hat{A}$ and $\hat{B}$ are homogeneous and coprime.

Introduce a new indeterminate $u$. Let $\sigma := (\sigma_0, \dots, \sigma_n) \in k^{n+1}$, and let

$$\tilde{F}_{\sigma}(x_0, \dots, x_n, u) := \hat{F}(x_0 u + \sigma_0, \dots, x_n u + \sigma_n),$$

$$\tilde{A}_{\sigma}(x_0, \dots, x_n, u) := \hat{A}(x_0 u + \sigma_0, \dots, x_n u + \sigma_n) := \tilde{A}_0 + \tilde{A}_1 u + \dots + \tilde{A}_{d_A} u^{d_A},$$

$$\tilde{B}_{\sigma}(x_0, \dots, x_n, u) := \hat{B}(x_0 u + \sigma_0, \dots, x_n u + \sigma_n) := \tilde{B}_0 + \tilde{B}_1 u + \dots + \tilde{B}_{d_B} u^{d_B},$$

where for $i = 0, \dots, d_A$ we have $\tilde{A}_i \in k[x_0, \dots, x_n]$ and for $i = 0, \dots, d_B$ we have $\tilde{B}_i \in k[x_0, \dots, x_n]$. If $\tilde{B}_0 \neq 0$, then $\tilde{A}_{\sigma}, \tilde{B}_{\sigma}$ are coprime by [124, Lemma 1]. In the following, assume $\tilde{B}_0 \neq 0$.

Let $\omega := (\omega_0, \dots, \omega_n) \in (k \setminus \{0\})^{n+1}$ be an admissible ratio (as in Section 2.4). For $i \in \mathbb{N}_0$, introduce notation:

$$\tilde{F}_{\sigma}^{[i]}(u) := \tilde{F}_{\sigma}(\omega_0^i, \dots, \omega_n^i, u),$$

$$\tilde{A}_{\sigma}^{[i]}(u) := \tilde{A}_{\sigma}(\omega_0^i, \dots, \omega_n^i, u) := \tilde{A}_0^{[i]} + \tilde{A}_1^{[i]} u + \dots + \tilde{A}_{d_A}^{[i]} u^{d_A},$$

$$\tilde{B}_{\sigma}^{[i]}(u) := \tilde{B}_{\sigma}(\omega_0^i, \dots, \omega_n^i, u) := \tilde{B}_0^{[i]} + \tilde{B}_1^{[i]} u + \dots + \tilde{B}_{d_B}^{[i]} u^{d_B}.$$

Assume $\tilde{A}_{\sigma}^{[i]}(u)$ and $\tilde{B}_{\sigma}^{[i]}(u)$ are coprime. Let $D := d_A+d_B+2$. Let $u_1, \dots, u_D \in k \setminus \{0\}$ be pairwise different. If $d_A$ and $d_B$ are known, then it is possible to obtain $\tilde{F}_{\sigma}^{[i]}(u)$ from the evaluations $\tilde{F}_{\sigma}^{[i]}(u_1), \dots, \tilde{F}_{\sigma}^{[i]}(u_D)$ using univariate rational interpolation (as in Section 2.3), assuming none of the evaluations fail. Then, for a sufficiently large $T$, it is possible to compute $\tilde{A}_{d_A}/\tilde{B}_0$ from $\tilde{A}_{d_A}^{[0]}/\tilde{B}_0, \dots, \tilde{A}_{d_A}^{[2T-1]}/\tilde{B}_0$ using sparse multivariate polynomial interpolation (as in Section 2.4). Division by $\tilde{B}_0$ ensures that the trailing coefficient of the denominator of $\tilde{F}_{\sigma}^{[i]}(u)$ is normalized to one across all $i = 0, \dots, 2T-1$. Similarly, it is possible to compute $\tilde{B}_{d_B}/\tilde{B}_0$.

Finally, note that $\tilde{A}_{d_A}, \tilde{B}_{d_B}$ coincide with $\hat{A}, \hat{B}$, since $\hat{A}$ and $\hat{B}$ are homogeneous. Then, by setting $x_0 := 1$, we obtain the sought $A = \hat{A}(1, x_1, \dots, x_n)$ and $B = \hat{B}(1, x_1, \dots, x_n)$.
**Algorithm 2** Sparse interpolation of a rational function given as a blackbox, [124, Section 3.3.]

**Input:** a blackbox BB(a) = F(a) or BB(a) returns FAIL; a subset S ⊆ k \ {0}; integers δA, δB; an admissible ratio ω ∈ (k \ {0})n+1.

**Output:** P, Q ∈ k[x] with P/Q = F and gcd(P, Q) = 1, or FAIL.

(Step 1) Let D := δA + δB + 2.

(Step 2) Choose the coordinates of σ := (σ₀, ..., σn) independently and uniformly from S ⊆ k \ {0}; choose pairwise different u₁, ..., uD ∈ k \ {0}.

(Step 3) For T = 1, 2, 4, ... do

(a) For i = 0, ..., 2T − 1 do

(i) Compute ˜Fσ [i](u₁), ..., ˜Fσ [i](uD) by evaluating the blackbox BB at the appropriate points; if any of the evaluations fail, **return** FAIL.

(ii) Compute ˜Fσ [i](u) by applying univariate rational interpolation to ˜Fσ [i](u₁), ..., ˜Fσ [i](uD) and (Dₐ, Dₐ) := (δA, δB) (as in Section 2.3); if the interpolation fails, **return** FAIL.

(b) Let P ∈ k[x₀, ..., xₙ] (resp., Q) be the output of sparse interpolation (Algorithm 1) applied to ˜A[0]dA / ˜B₀, ..., ˜A[2T − 1]dA / ˜B₀ (resp., ˜B[0]dB / ˜B₀, ..., ˜B[2T − 1]dB / ˜B₀).

(c) If (δA, δB) = (deg P, deg Q) and P(σ) = Q(σ) ˜F (σ), where ˜F (σ) is computed by evaluating the blackbox BB, then **return** P(1, x₁, ..., xₙ) and Q(1, x₁, ..., xₙ).

(d) If T > (n+δA+δB) / 2, then **return** FAIL.

**Remark 2.7.** One difference between the presentation here and in [124, Section 3.3.] is that here we allow for the input (δA, δB) to be different from the actual (dA, dB) (and the algorithm will either fail or produce an incorrect result in this case). The condition in (Step 3)(d) is to make sure that the algorithm terminates when δA ≠ dA or δB ≠ dB.

**Proposition 2.8.** Let g ∈ k[x] be a nonzero polynomial such that for all a ∈ kⁿ where BB(a) returns FAIL we have g(a) = 0. Let h := deg g. Let s := max(sA, sB) and D := dA + dB + 2.

Algorithm 2 terminates for any inputs. If δA = dA and δB = dB, then Algorithm 2 does not fail and returns a correct result with probability at least

$$1 - \frac{D(\lceil \log_2 s \rceil + 8s + 4sh + 8sD)}{|S|}.$$

*Proof.* Correctness is proven in [124, Theorem 2]. For the termination claim, note that the algorithm terminates for any inputs in (Step 3)(d) for large enough T, in the case it did not terminate prior to that.

We treat probability in the case δA = dA and δB = dB. Following Remark 2.6, we shall replace ˜Fσ (x₀, ..., xₙ, u) with ˜Fσ (γ₀x₀, ..., γₙxₙ, u), where the coordinates of γ := (γ₀, ..., γₙ) are chosen independently and uniformly from S ⊆ k \ {0}. The algorithm terminates no later than the first iteration when T ≥ s. Therefore, upon termination, we have T ≤ 2s. Algorithm 2 returns FAIL only in one of the following cases:

• We have ˜B₀ = 0. This is equivalent to ˜B(σ₀, ..., σₙ) = 0. This can be checked when performing rational function interpolation in (Step 3)(a)(ii), assuming the interpolation does not fail. Note deg ˜B ≤ D. Therefore, the probability that ˜B(σ₀, ..., σₙ) = 0 does not exceed D/|S| by Lemma 2.1.

• For any of i = 0, ..., 2T − 1, any of the evaluations ˜Fσ [i](u₁), ..., ˜Fσ [i](uD) fail. For every i = 0, ..., 2T − 1 and j = 1, ..., D, evaluating ˜Fσ [i](uⱼ) is equivalent to evaluating BB(ξ [i,j]) with

$$\xi^{[i,j]} := \left( \frac{\gamma_1 \omega_1^i u_j + \sigma_1}{\gamma_0 \omega_0^i u_j + \sigma_0}, \dots, \frac{\gamma_n \omega_n^i u_j + \sigma_n}{\gamma_0 \omega_0^i u_j + \sigma_0} \right) \in \mathbb{K}^n.$$

Failure occurs when either of the two happens: a division by zero when constructing ξ [i,j], or the blackbox returning failure. The probability of the former is at most 4sD/|S|, which can be seen by applying the Schwartz-Zippel lemma (Lemma 2.1) to p ∈ k[σ] with deg p ≤ 4sD and the evaluation point σ₀, where

$$p := \prod_{\substack{i=0,\dots,2T-1 \\ j=1,\dots,D}} (\gamma_0 \omega_0^i u_j + \sigma).$$
For the latter, recall that for all $\mathbf{a} \in k^n$ where the blackbox fails we have $g(\mathbf{a}) = 0$ with $h = \deg g$. We apply Corollary 2.2 to $g$ and $\xi^{[i,j]}$, which is justified since none of $\omega_1, \dots, \omega_n$ and $u_1, \dots, u_D$ are zero. Whence, and because of $2T \le 4s$, the probability of failure does not exceed $4shD/|S|$.

* For any of $i = 0, \dots, 2T - 1$, rational univariate interpolation of $\tilde{F}_{\sigma}^{[i]}(u)$ fails. This can be detected by comparing the degrees of the results against $(\delta_A, \delta_B)$. As shown in [124, Section 3.3.], there exists a nonzero polynomial $R_{\sigma} \in k[x_0, \dots, x_n]$ with $\deg R_{\sigma} \le 2D^2$ such that for every $i = 0, \dots, 2T - 1$ rational univariate interpolation of $\tilde{F}_{\sigma}^{[i]}(u)$ fails whenever $R_{\sigma}(\gamma\omega^i) = 0$. Since $2T \le 4s$, the probability of a failure does not exceed $8sD^2/|S|$ by Corollary 2.2.

Provided Algorithm 2 does not return FAIL, result may be incorrect. This happens when the check in (Step 3)(c) passes for incorrect $P, Q$. Note $P, Q$ are correct if and only if $\deg P = \deg A$ and $\deg Q = \deg B$ and $P = Q\hat{F}$. Therefore, it suffices to estimate the probability that $P\hat{B} - Q\hat{A} \ne 0$ vanishes at $\sigma$. We apply the Schwartz-Zippel lemma (Lemma 2.1), which is well-defined because the evaluation $\hat{F}(\sigma)$ does not fail since $\hat{B}_0 \ne 0$. There are at most $\lceil \log_2 s \rceil$ iterations where the Schwartz-Zippel lemma is applied, thus the probability of producing an incorrect result is bounded by $D[\log_2 s]/|S|$. Putting everything together, we arrive at the desired probability. $\square$

## 2.6 The degree of a blackbox function

Let $F = A/B$ be a rational function in $k(\mathbf{x})$ with $A, B \in k[\mathbf{x}]$ and $\gcd(A, B) = 1$, with $\mathbf{x} = (x_1, \dots, x_n)$. As before, assume we have access to a blackbox that evaluates $F$ at a point in $k^n$. Let $d_A := \deg A$ and $d_B := \deg B$. We describe an algorithm for computing $d_A$ and $d_B$.

In the notation of the previous section, let $\sigma := (\sigma_0, \dots, \sigma_n) \in k^{n+1}$ be arbitrary and let

$$R_{\sigma} := \mathrm{Res}_{u}(\tilde{A}_{\sigma}(x_0, \dots, x_n, u), \tilde{B}_{\sigma}(x_0, \dots, x_n, u)),$$

where $\mathrm{Res}_u$ denotes the resultant with respect to the variable $u$. Note this is the same $R_{\sigma}$ that appears in the proof of Proposition 2.8. By [124, Theorem 2], $R_{\sigma}$ is not identically zero. By definition, for every $(a_0, \dots, a_n) \in k^{n+1}$, the degrees of the numerator and denominator of $\tilde{F}_{\sigma}(a_0, \dots, a_n, u)$ coincide with that of $F$ if and only if $R_{\sigma}(a_0, \dots, a_n) \ne 0$.

For convenience, introduce notation

$$\tilde{f}_{\sigma}(u) := \tilde{F}_{\sigma}(1, \dots, 1, u).$$

When $R_{\sigma}(1, \dots, 1) \ne 0$, the degrees of the numerator and denominator of $\tilde{f}_{\sigma}$ and $\tilde{F}_{\sigma}$ coincide, and are equal to $d_A$ and $d_B$, respectively. Let $T \ge \max(d_A, d_B)$ and let $u_1, \dots, u_{2T+2} \in k \setminus \{0\}$ be pairwise different. Then, we can recover $\tilde{f}_{\sigma}$ by applying univariate rational interpolation to the evaluations $\tilde{f}_{\sigma}(u_1), \dots, \tilde{f}_{\sigma}(u_{2T+2})$ and the degree bounds $D_A := D_B := T$ (as in Section 2.3). This leads to the following algorithm.

**Algorithm 3** The degree of a rational function given as a blackbox

**Input:** a blackbox BB, such that for every $\mathbf{a} \in k^n$ either $BB(\mathbf{a}) = F(\mathbf{a})$ or $BB(\mathbf{a})$ returns FAIL; an integer $d$; a subset $S \subseteq k \setminus \{0\}$.

**Output:** if $d_A + d_B \le d$, the integers $d_A, d_B$ or FAIL; if $d_A + d_B > d$, STOPPED or FAIL.

(Step 1) Choose the coordinates of $\sigma := (\sigma_0, \dots, \sigma_n)$ independently and uniformly from $S \subseteq k \setminus \{0\}$.

(Step 2) Choose $u_0$ uniformly from $S \subseteq k \setminus \{0\}$. If the evaluation $\tilde{f}_{\sigma}(u_0)$ fails (computed by evaluating the blackbox BB), return FAIL.

(Step 3) For $T = 0, 1, 2, \dots$ do

(a) If $T > d$, then return STOPPED.

(b) Compute $\tilde{f}_{\sigma}(u_1), \dots, \tilde{f}_{\sigma}(u_{2T+2})$ for pairwise different $u_1, \dots, u_{2T+2} \in k \setminus \{0\}$ by evaluating the blackbox BB at the appropriate points; if any of the evaluations fail, return FAIL.

(c) Compute $p/q \in k(u)$ with $\gcd(p, q) = 1$ by applying univariate rational interpolation to $\tilde{f}_{\sigma}(u_1), \dots, \tilde{f}_{\sigma}(u_{2T+2})$ with the degree bounds $D_A := D_B := T$ (as in Section 2.3).

(d) If $p(u_0) = q(u_0)\tilde{f}_{\sigma}(u_0)$, where $\tilde{f}_{\sigma}(u_0)$ is computed by evaluating the blackbox BB, then return $d_p, d_q$ if $d_p + d_q \le d$ and STOPPED, otherwise.
**Proposition 2.9.** Let $g \in k[x]$ be a nonzero polynomial such that for all $\mathbf{a} \in k^n$ where $\mathrm{BB}(\mathbf{a})$ returns FAIL we have $g(\mathbf{a}) = 0$. Let $h := \deg g$. Let $D := d_A + d_B + 2$.

Algorithm 3 terminates for any inputs. The probability that Algorithm 3 does not fail and returns a correct result is at least

$$1 - \frac{D(2 + 3h + 3D)}{|S|}.$$

*Proof.* Termination is clear due to (Step 3)(a).

We prove correctness and probability. Following Remark 2.6, we replace $\tilde{F}_{\sigma}(x_0, \dots, x_n, u)$ with $\tilde{F}_{\sigma}(\gamma_0 x_0, \dots, \gamma_n x_n, u)$, where the coordinates of $\gamma := (\gamma_0, \dots, \gamma_n)$ are chosen independently and uniformly from $S \subseteq k \setminus \{0\}$. Note that at every iteration we have $T + 1 \le D$. In the case $d_A + d_B > d$, this is clear; in the case $d_A + d_B \le d$, this is because the algorithm terminates no later than the first iteration when $T = \max(d_A, d_B)$.

The algorithm returns FAIL only in one of the following cases:

* Any of the evaluations $\tilde{f}_{\sigma}(u_1), \dots, \tilde{f}_{\sigma}(u_{2T+2})$ fail. Similarly to the proof of Proposition 2.8, failure occurs either because of division by zero when constructing the evaluation point or when the evaluation of the blackbox returns fail.

The probability of these does not exceed $(2D + 2hD)/|S|$ by Corollary 2.2.

* The evaluation $\tilde{f}_{\sigma}(u_0)$ fails. By the Schwartz-Zippel lemma (Lemma 2.1), the probability of this does not exceed $h/|S|$.

In the following, assume the algorithm did not fail. We treat the case $d_A + d_B > d$. Incorrect result is returned only if the check in (Step 3)(d) ever passes. If $R_{\sigma}(1, \dots, 1) \ne 0$, then at every iteration $p - q\tilde{f}_{\sigma} \ne 0$. We appeal to Remark 2.6, so it suffices to estimate the probability that either $R_{\sigma}(\gamma_0, \dots, \gamma_n)$ vanishes or that at any iteration $p - q\tilde{f}_{\sigma}$ vanishes at $u_0$. The former is at most $2D^2/|S|$ by the Schwartz-Zippel lemma (Lemma 2.1), since $\deg R_{\sigma} \le 2D^2$. The latter is at most $D^2/|S|$, since we apply the Schwartz-Zippel lemma (Lemma 2.1) at most $D$ times with the same $u_0$, each time to a polynomial of degree at most $D$.

When the algorithm does not fail and does not return an incorrect result, the output is correct. Thus, we obtain the desired probability of correctness in the case $d_A + d_B > d$. The case $d_A + d_B \le d$ is similar: the output in (Step 3)(d) is guaranteed to be correct at the first iteration when $T = \max(d_A, d_B)$ in the case $R_{\sigma}(1, \dots, 1) \ne 0$ and assuming the algorithm did not fail and did not terminate with incorrect result. $\square$

## 2.7 OMS Ideals

Let $k(\mathbf{x})$ be the field of rational functions in the indeterminates $\mathbf{x} := (x_1, \dots, x_n)$. We will be interested in subfields $k \subset \mathcal{E} \subseteq k(\mathbf{x})$. By [84, Ch. VIII, Ex. 4] every such subfield is finitely generated over $k$, that is, $\mathcal{E} = k(\mathbf{g})$ for some $\mathbf{g} = (g_1, \dots, g_m)$ with $g_1, \dots, g_m \in k(\mathbf{x})$.

One of the key tools in this paper are certain ideals associated to subfields $k \subset \mathcal{E} \subseteq k(\mathbf{x})$. These ideals were introduced by Ollivier in [104, Ch. III, §2] and Müller-Quade and Steinwandt in [101]. Thus, we will call them OMS ideals.

**Definition 2.10 (OMS ideal).** Let $\mathcal{E} \supset k$ be a subfield of $k(\mathbf{x})$ represented as $\mathcal{E} = k(g_1, \dots, g_m)$, where $g_i = \frac{p_i(\mathbf{x})}{q_i(\mathbf{x})}$ for $i = 1, \dots, m$. We denote by $Q$ the $\mathrm{lcm}$ of $q_1, \dots, q_m$. Introduce new indeterminates $\mathbf{y} := (y_1, \dots, y_n)$. The following ideal in $k(\mathbf{x})[\mathbf{y}]$ will be called the OMS ideal of $\mathcal{E}$:

$$\mathrm{OMS}_{\mathcal{E}} := \langle p_i(\mathbf{y})q_i(\mathbf{x}) - q_i(\mathbf{y})p_i(\mathbf{x}) \mid 1 \le i \le m \rangle : (Q(\mathbf{y}))^{\infty}. \quad (3)$$

**Remark 2.11 (Dependence on the choice of generators).** Formally, the definition (3) depends not only on $\mathcal{E}$ but also on the chosen set of generators. In fact, [101, Lemma 1.2] shows that $\mathrm{OMS}_{\mathcal{E}}$ is generated by $\langle y_1 - x_1, \dots, y_n - x_n \rangle \cap \mathcal{E}[\mathbf{y}]$ and, thus, does not depend on the choice of generators.

This also implies that, if $\mathrm{char} k = 0$, then $\mathrm{OMS}_{\mathcal{E}}$ is radical. It may be not the case in positive characteristic: $\mathrm{OMS}_{\mathbb{F}_p(x^p)}$ is not radical.

OMS ideals are useful for characterizing the elements of rational function fields algorithmically by associating rational functions with polynomials of special form:

**Lemma 2.12 ([101, Lemmas 1.2 and 1.6.], [104, p. 51, Theorem 3]).** Let $\mathcal{E} \supset k$ be a subfield of $k(\mathbf{x})$. Then:

1. for every polynomials $p, q \in k[x]$, $\mathfrak{p} \in \mathcal{E}$ if and only if $p(y)q(x) - q(y)p(x) \in \mathrm{OMS}_{\mathcal{E}}$.
2. Let $f_1, \dots, f_s$ be the reduced Gröbner basis of OMS$_\mathcal{E}$ with respect to any ordering. Then the coefficients of $f_1, \dots, f_s$ generate $\mathcal{E}$ over $k$.

The second property from the lemma provides, in particular, a way to compute a new (often, simpler) set of generators for $\mathcal{E}$ starting from some generators.

**Example 2.13.** Let $\mathcal{E} = \mathbb{Q}(x_1^2 + x_2^2, x_1^3 + x_2^3, x_1^4 + x_2^4) \subset \mathbb{Q}(x_1, x_2)$. We have $Q(y_1, y_2) = 1$ and

$$ \mathrm{OMS}_{\mathcal{E}} = \left\{ \begin{array}{l} y_1^2 + y_2^2 - (x_1^2 + x_2^2), \\ y_1^3 + y_2^3 - (x_1^3 + x_2^3), \\ y_1^4 + y_2^4 - (x_1^4 + x_2^4) \end{array} \right\} \subset \mathbb{Q}(x_1, x_2)[y_1, y_2]. $$

The reduced Gröbner basis of OMS$_\mathcal{E}$ with respect to the lexicographic ordering on $y_1, y_2$ with $y_2 < y_1$ is

$$ y_1 + y_2 - (x_1 + x_2), \quad y_2^2 - (x_1 + x_2) + x_1 x_2. $$

By Lemma 2.12, the underlined coefficients generate $\mathcal{E}$, so $\mathcal{E} = \mathbb{Q}(x_1 + x_2, x_1 x_2)$. Thus, the computation of the Gröbner basis of the OMS ideal allowed us to recognize that $\mathcal{E}$ is the field of symmetric functions and provided us with the most natural set of generators.

Since the definition (3) involves saturation, practical computation with OMS$_\mathcal{E}$ requires passing through an elimination ordering (if the saturation is done via the Rabinowitsch trick) or using dedicated algorithms for computing Gröbner bases of saturations. It would naturally simplify the computations if we could employ the Rabinowitsch trick with respect to a new variable $t$ without performing elimination.

**Notation 2.14 (eOMS).** In the notation of Definition 2.10, denote

$$ \mathrm{eOMS}_g := \langle p_1(\mathbf{y})q_1(\mathbf{x}) - q_1(\mathbf{y})p_1(\mathbf{x}), \dots, p_m(\mathbf{y})q_m(\mathbf{x}) - q_m(\mathbf{y})p_m(\mathbf{x}), tQ(\mathbf{y}) - 1 \rangle \subset k(\mathbf{x})[\mathbf{y}, t]. \quad (4) $$

Obviously, the first statement from Lemma 2.12 holds if we replace OMS$_\mathcal{E}$ with eOMS$_g$. As we prove in the following lemma, the same is true for the second statement as well.

**Lemma 2.15.** Let $\mathcal{E} \supset k$ be a subfield of $k(\mathbf{x})$ represented by $\mathcal{E} = k(\mathbf{g})$. Let $f_1, \dots, f_s$ be the reduced Gröbner basis of eOMS$_g$ with respect to any ordering (see (4)). Then the coefficients of $f_1, \dots, f_s$ generate $\mathcal{E}$ over $k$.

*Proof.* Let $\mathcal{F}$ be the field generated over $k$ by the coefficients of $f_1, \dots, f_s$. Since the original ideal generators from (4) have coefficients in $\mathcal{E}$ and Gröbner basis computation does not extend the ground field, we conclude that $\mathcal{F} \subset \mathcal{E}$. Assume that $\mathcal{F} \neq \mathcal{E}$. Without loss of generality, $g_1 \notin \mathcal{F}$. We choose a basis $\{e_i\}_{i \in \Lambda}$ of $\mathcal{E}$ over $\mathcal{F}$, where $\Lambda$ is either $\{1, 2, \dots, \dim_{\mathcal{F}} \mathcal{E}\}$ if $\dim_{\mathcal{F}} \mathcal{E} < \infty$ or $\mathbb{N}$ otherwise. The basis can be chosen so that $e_1 = 1$. We write $g_1$ in this basis as $\sum_{i \in \Lambda} c_i e_i$, where $c_i \in \mathcal{F}$ for every $i \in \Lambda$. Since $g_1 \notin \mathcal{F}$,

there exists $i_0 \in \Lambda \setminus \{1\}$ such that $c_{i_0} \neq 0$. We write $g_1 = \frac{p_1(\mathbf{x})}{q_1(\mathbf{x})}$ for $p_1, q_1 \in k[\mathbf{x}]$. Let $r_1$ and $r_2$ be the normal forms of $p_1(\mathbf{y}) - c_1(\mathbf{x})q_1(\mathbf{y})$ and $q_1(\mathbf{y})$ with respect to $f_1, \dots, f_s$. Since reductions do not extend the ground field, $r_1, r_2 \in \mathcal{F}[\mathbf{y}, t]$. Using the expansion for $g_1$ with respect to the basis, $p_1(\mathbf{y}) - g_1(\mathbf{x})q_1(\mathbf{y})$ can be written as

$$ (p_1(\mathbf{y}) - c_1(\mathbf{x})q_1(\mathbf{y}))e_1 - \sum_{i \in \Lambda \setminus \{1\}} c_i(\mathbf{x})q_1(\mathbf{y})e_i. $$

Then the normal form of $p_1(\mathbf{y}) - g_1(\mathbf{x})q_1(\mathbf{y})$ with respect to $f_1, \dots, f_s$ is equal to

$$ r_1 e_1 - \sum_{i \in \Lambda \setminus \{1\}} c_i r_2 e_i. $$

This expression is equal to zero, so we have $c_i r_2 e_i = 0$ for every $i \in \Lambda \setminus \{1\}$. In particular, $c_{i_0} r_2 e_{i_0} = 0$, so $r_2 = 0$. Then $q_1(\mathbf{y}) \in \mathrm{eOMS}_g \implies Q(\mathbf{y}) \in \mathrm{eOMS}_g \implies 1 \in \mathrm{eOMS}_g$. This contradicts the fact that eOMS$_g$ has a solution $\mathbf{y} = \mathbf{x}$. The contradiction implies that $\mathcal{F} = \mathcal{E}$ as desired. $\square$

**Example 2.16 (OMS may bear excess).** While coefficients of the OMS ideal turn out to be very useful in finding simpler generating sets for subfields of $k(\mathbf{g})$, not all the coefficients are equally useful. Fix a positive integer $\ell$ and consider a subfield in $k(x_1, x_2)$ generated by $g_1 = x_1 + x_2$ and $g_2 = \frac{x_1^\ell}{x_2^\ell}$. Then one can show that the degrevlex Gröbner basis of OMS$_{k(g_1, g_2)}$ will consist of $y_1 + y_2 - (x_1 + x_2)$ (which gives one generator, $x_1 + x_2$) and the following polynomial:

$$ y_2^\ell - \underbrace{\frac{\ell x_2^\ell (x_1 + x_2)}{x_2^\ell + (-1)^{\ell+1} x_1^\ell}}_{\text{second generator}} y_2^{\ell-1} + \sum_{k=0}^{\ell-2} (-1)^{\ell-k} \binom{\ell}{k} \underbrace{\frac{(x_1 + x_2)^{\ell-k-1} b^\ell}{x_2^\ell + (-1)^{\ell+1} x_1^\ell}}_{\text{useless coefficients}} y_2^k. $$
Not only the majority of the coefficients are not needed to generate the field but also the degree of the numerator of the extra coefficients may reach $2\ell - 1$, which is almost twice larger than that of the original generators. One observes similar phenomenon when replacing $\text{OMS}_{k(g_1,g_2)}$ with $\text{eOMS}_{g_1,g_2}$.

**Example 2.17 (Dura lex, sed lex).** The generating set coming from the coefficients of the OMS ideal may depend significantly on the ordering used for the Gröbner basis computation. Consider a tuple $g := (x_1(x_2-x_3), x_2+x_3, x_2x_3)$ of elements of $k(x_1, x_2, x_3)$. Computing the Gröbner basis of $\text{OMS}_{k(g)}$ with respect to either the lexicographic or the degrevlex ordering with $y_1 > y_2 > y_3$ yields

$$y_3^2 - (x_2+x_3)y_3 + x_2x_3, \quad y_2+y_3 - (x_2+x_3), \quad y_1+2\frac{x_1}{x_2-x_3}y_3 - \frac{x_1(x_2+x_3)}{x_2-x_3}.$$

This basis gives back the original generators $x_2+x_3$ and $x_2x_3$ and provides additional $\frac{x_1}{x_2-x_3}$ and $\frac{x_1(x_2+x_3)}{x_2-x_3}$, which do not look simpler than the original $x_1(x_2-x_3)$. On the other hand, the Gröbner basis of $\text{OMS}_{k(g)}$ with respect to either the lexicographic or the degrevlex ordering with $y_3 > y_2 > y_1$ is equal to

$$y_1^2 - x_1^2, \quad y_2 + \frac{x_3 - x_2}{2x_1}y_1 - \frac{x_2 + x_3}{2}, \quad y_3 + \frac{x_2 - x_3}{2x_1}y_1 - \frac{x_2 + x_3}{2}.$$

Here, a simple generator $x_2x_3$ has disappeared but another simple generator, namely $x_1^2$, is discovered.

## 2.8 Randomized subfield membership testing

**Lemma 2.12** provides a constructive way to check if a given function $f \in k(\mathbf{x})$ belongs to a subfield $k(g_1, \dots, g_m)$. However, this approach would involve computing Gröbner basis with coefficients in the field of rational functions $k(\mathbf{x})$. This may be quite heavy and one of the leitmotifs of the present paper is to avoid such computations by working with specializations of the corresponding ideals.

**Notation 2.18 (Specialization of eOMS).** In the notation of Definition 2.10, consider $\mathbf{a} \in k^n$. We define $\text{eOMS}_{\mathbf{a}}(\mathbf{a})$ by:

$$\text{eOMS}_{\mathbf{a}}(\mathbf{a}) := \langle p_1(\mathbf{y})q_1(\mathbf{a}) - q_1(\mathbf{y})p_1(\mathbf{a}), \dots, p_m(\mathbf{y})q_m(\mathbf{a}) - q_m(\mathbf{y})p_m(\mathbf{a}), tQ(\mathbf{y}) - 1 \rangle \subseteq k[\mathbf{y}, t].$$

For membership testing, we will rely on the following fact.

**Theorem 2.19 ([37, Theorem 3.3]).** Let $k$ be a field of characteristic zero. Let $p, p_1, \dots, p_m, q \in k[\mathbf{x}]$, where $\mathbf{x} = (x_1, \dots, x_n)$. Let $0 < \varepsilon < 1$ be a real number. We define

$$d := \max(\deg q + 1, \deg p, \deg p_1, \dots, \deg p_m) \quad \text{and} \quad M := \frac{6d^{n+3}}{\varepsilon}.$$

Let $\mathbf{a} := (a_1, \dots, a_n)$ be integers sampled uniformly and independently from $[0, M]$. Consider the field $\mathcal{E} := k\left(\frac{p_1(\mathbf{x})}{q(\mathbf{x})}, \dots, \frac{p_m(\mathbf{x})}{q(\mathbf{x})}\right)$. Then,

1. if $p(\mathbf{y})q(\mathbf{a}) - p(\mathbf{a})q(\mathbf{y}) \in \text{eOMS}_{\mathbf{a}}(\mathbf{a})$, then $\frac{p(\mathbf{x})}{q(\mathbf{x})} \in \mathcal{E}$ with probability at least $1 - \varepsilon$;

2. if $p(\mathbf{y})q(\mathbf{a}) - p(\mathbf{a})q(\mathbf{y}) \notin \text{eOMS}_{\mathbf{a}}(\mathbf{a})$, then $\frac{p(\mathbf{x})}{q(\mathbf{x})} \notin \mathcal{E}$ with probability at least $1 - \varepsilon$.

This theorem yields a randomized Monte-Carlo algorithm for checking field membership requiring Gröbner basis computation only over $k$. We can go a bit further and use the approach from [68] to ensure that the ideal for which the Gröbner basis computation is performed is of dimension zero. This is based on the following lemma implicitly contained in [68, Proof of Theorem 1] (also reminiscent to the use of cross-section in rational invariant computation [64, Section 3]).

**Lemma 2.20.** Consider a subfield $\mathcal{E} = k(g_1, \dots, g_m) \subseteq k(\mathbf{x})$, and let $z_1, \dots, z_s$ be a transcendence basis of $k(\mathbf{x})$ over $\mathcal{E}$. Assume that $g \in k(\mathbf{x})$ is algebraic over $\mathcal{E}$. Then $g \in \mathcal{E} \iff g \in \mathcal{E}(z_1, \dots, z_s)$.

*Proof.* The implication $g \in \mathcal{E} \implies g \in \mathcal{E}(z_1, \dots, z_s)$ is trivial, let us prove the converse. Assume that $g \in \mathcal{E}(z_1, \dots, z_s)$ and $g$ is algebraic over $\mathcal{E}$. Let $g_1, \dots, g_\ell$ be all the roots of the minimal polynomial for $g$ over $\mathbb{K}(\mathbf{x})$ such that $g_1 = g$. We write $g = G(z_1, \dots, z_s)$, where $G$ is a rational function from $\mathcal{E}(Z_1, \dots, Z_s)$ with $Z_1, \dots, Z_s$ being new indeterminates. Assume that $G$ is nonconstant and, without loss of generality, that it depends nontrivially on $Z_1$. Then $G(z_1 + X, z_2, \dots, z_s)$ is a nonconstant rational function in $X$. Since $\bar{k}$ is infinite, there exists $c \in \bar{k}$ such that $G(z_1 + c, z_2, \dots, z_s) \neq g_i$ for every $1 \le i \le \ell$. Since $z_1, \dots, z_s$ are transcendental over $\mathcal{E}$, there exists [95, Proposition 2.4(a)] an automorphism $\alpha: \bar{k}(\mathbf{x}) \to \bar{k}(\mathbf{x})$ such that $\alpha|_{\mathcal{E}}$ is the identity, $\alpha(z_1) = z_1 + c$ (where $c$ is viewed as an element of $\bar{k} \subset \bar{k}(\mathbf{x})$), and $\alpha(z_i) = z_i$ for every $i = 2, \dots, s$. Then $\alpha(g)$ satisfies the same equation over $\mathcal{E}$ as $g$, so it must belong to the set $\{g_1, \dots, g_\ell\}$. This contradicts the choice of $c$. Therefore, $G \in \mathcal{E}$, so $g \in \mathcal{E}$. $\square$
Combining Theorem 2.19 and Lemma 2.20, we arrive at the following algorithm.

**Algorithm 4 Randomized subfield membership check**

**Input:** polynomials $p, p_1, \dots, p_m, q \in k[x]$, where $\text{char } k = 0$, and a real number $0 < \varepsilon < 1$;

**Output:** True if $\frac{p}{q} \in k(\frac{p_1}{q}, \dots, \frac{p_m}{q})$, and False otherwise.

(Step 1) Set $d := \max(\deg q + 1, \deg p, \deg p_1, \dots, \deg p_m)$.

(Step 2) Compute the Jacobian $J(\mathbf{x})$ of $\frac{p_1}{q}, \dots, \frac{p_m}{q}$ (we assume that the gradients are the rows of $J(\mathbf{x})$) and the gradient $G(\mathbf{x})$ of $\frac{p}{q}$.

(Step 3) Sample each coordinate of $\mathbf{a} \in \mathbb{Z}^n$ uniformly at random from $[0, \frac{2(2 \min(m,n)+1)(d-1)}{\varepsilon}] \cap \mathbb{Z}$.

(Step 4) If $G(\mathbf{a})$ does not belong to the rowspan of $J(\mathbf{a})$, return False.

(Step 5) Take any set of columns of $J(\mathbf{a})$ forming a basis of the column space of the matrix. Let $S \subset \{1, \dots, n\}$ be the indices these columns. Set $\bar{S} := \{1, \dots, n\} \setminus S$.

(Step 6) Sample $\mathbf{b} \in \mathbb{Z}^n$ with each coordinate drawn uniformly at random from $[0, 12d^{n+3}/\varepsilon] \cap \mathbb{Z}$.

(Step 7) Using Gröbner basis, check if $p(\mathbf{x})q(\mathbf{b}) - p(\mathbf{b})q(\mathbf{x})$ belongs to

$$e\text{OMS}_{\mathbf{p}/q}(\mathbf{b}) + \langle x_j - b_j \mid j \in \bar{S} \rangle.$$

(Step 8) If belongs, return True, otherwise return False.

**Proposition 2.21.** *The output of Algorithm 4 is correct with the probability at least $1 - \varepsilon$.*

*Proof.* Denote $g := \frac{p}{q}, g_1 := \frac{p_1}{q}, \dots, g_m := \frac{p_m}{q}$, and $\mathcal{E} := k(g_1, \dots, g_m)$. Let $\tilde{J}(\mathbf{x})$ be the matrix obtained by adding an extra row $G(\mathbf{x})$ to $J(\mathbf{x})$.

First, we consider the case when $g$ is transcendental over $\mathcal{E}$. By [38, Proposition 2.4], this implies $\text{rank } J(\mathbf{x}) < \text{rank } \tilde{J}(\mathbf{x})$ and, in particular, any maximal nonsingular minor of $\tilde{J}(\mathbf{x})$ contains the row corresponding to $G$. Let $M(\mathbf{x})$ be any maximal nonsingular minor of $\tilde{J}(\mathbf{x})$. Assume that $q(\mathbf{a}) \neq 0$ (so that $\tilde{J}(\mathbf{a})$ is well-defined) and $\det M(\mathbf{a}) \neq 0$. Then

$$\text{rank } J(\mathbf{a}) \le \text{rank } J(\mathbf{x}) < \text{rank } \tilde{J}(\mathbf{x}) = \text{rank } \tilde{J}(\mathbf{a}),$$

so the algorithm will return False at (Step 4) which will be the correct answer. The degree of the numerator of $\det M(\mathbf{x})$ does not exceed $\min(m,n)(2d-2)$. Therefore, by the Schwartz-Zippel lemma (Lemma 2.1), $\tilde{J}(\mathbf{x})$ will be well-defined under $\mathbf{x} \to \mathbf{a}$ and the correct result will be returned with the probability at least

$$1 - \frac{d-1 + \min(m, n)(2d-2)}{2(2 \min(m, n) + 1)(d-1)/\varepsilon} = 1 - \frac{\varepsilon}{2} > 1 - \varepsilon.$$

Next, we consider the case when $g$ is algebraic over $\mathcal{E}$. In this case, $G(\mathbf{x})$ belongs to the rowspan of $J(\mathbf{x})$ by [38, Proposition 2.4]. Let $M(\mathbf{x})$ be a maximal nonsingular minor of $J(\mathbf{x})$. Assume that $q(\mathbf{a}) \neq 0$ and $\det M(\mathbf{a}) \neq 0$. Then $\text{rank } J(\mathbf{a}) = \text{rank } J(\mathbf{x}) = \text{rank } \tilde{J}(\mathbf{x}) \ge \text{rank } \tilde{J}(\mathbf{a})$. Therefore, the algorithm will not return on (Step 4).

Furthermore, since $\text{rank } J(\mathbf{a}) = \text{rank } J(\mathbf{x})$, every set of columns which is a basis of the columns space of $J(\mathbf{a})$ is a specialization of a column basis for $J(\mathbf{x})$. Thus, [38, Proposition 2.4] implies that $\{x_j \mid j \in \bar{S}\}$ is a transcendence basis of $k(\mathbf{x})$ over $\mathcal{E}$. Then, by Lemma 2.20, $g \in \mathcal{E}$ if and only if $g \in \mathcal{E}(x_j \mid j \in \bar{S})$. The latter containment will be assessed by (Step 7) with the error probability not exceeding $\varepsilon/2$ by Theorem 2.19 (where the saturation ideal is replaced by the one given by the Rabinowitsch trick, this does not affect the membership). Using the same calculation as in the first part of the proof, the probability of $q(\mathbf{x}) \det M(\mathbf{x})$ vanishing at $\mathbf{a}$ also does not exceed $\varepsilon/2$. Therefore, the returned result will be correct with probability at least $1 - \varepsilon$. $\square$

**Remark 2.22.** In practice, in (Step 5), we cache the set $\bar{S}$ for a field $\mathcal{E}$ to reuse in further membership tests.

**Remark 2.23.** In some parts of our algorithm, we will not need the explicit error bound, so we will perform computation modulo large prime (and sampling from the whole field) instead of $k = \mathbb{Q}$. Then, for every input polynomials and every error probability $\varepsilon$, there will be a large enough prime (e.g., exceeding all the integers arising during the computation over $\mathbb{Q}$) so that the error bound holds for the modular version of Algorithm 4 as well.
**Remark 2.24 (On saturation).** When computing Gröbner bases of specialized eOMS ideals using our implementation (Section 5.5), we found that factoring the saturation polynomial $Q(y)$ (as in [101, Lemma 1.3.]) does not provide a measurable speed-up on our examples.

**Remark 2.25.** Algorithm 4 yields a randomized algorithm for testing equality of subfields of $k(x)$. Let $\mathcal{E} = k(f_1, \dots, f_s)$ and $\mathcal{F} = k(g_1, \dots, g_r)$ be two such subfields, and let $0 < \varepsilon < 1$ be a real number. Then we apply Algorithm 4 to check whether

$$ (\forall 1 \le i \le s: f_i \in \mathcal{F}) \ \& \ (\forall 1 \le j \le r: g_j \in \mathcal{E}). $$

If each individual check is performed with the probability of correctness at least $1 - \frac{\varepsilon}{s+r}$, then the result is correct with probability at least $1 - \varepsilon$.

# 3 Building blocks of the algorithm

## 3.1 Basic minimization

We start with a method that follows straightforwardly from the field membership algorithm. Let $g \subset k(x)$. The following method computes a subset $h \subseteq g$ that generates the same field as $g$ over $k$ and is minimal with respect to inclusion.

**Algorithm 5 Computing a minimal subset**

**Input:** $g \subset k(x)$, where $\text{char } k = 0$, and a real number $0 < \varepsilon < 1$;

**Output:** $h \subseteq g$ such that $k(h) = k(g)$ and $h$ is minimal with respect to inclusion.

(Step 1) Set $h := (h_1, \dots, h_m) := g$.

(Step 2) For $i = 1, \dots, m$ do

(a) If $h_i \in k(h \setminus \{h_i\})$ with probability at least $1 - \frac{\varepsilon}{m}$ using Algorithm 4, then set $h := h \setminus \{h_i\}$.

(Step 3) Return $h$.

**Proposition 3.1.** *Algorithm 5 is correct with probability at least $1 - \varepsilon$.*

*Proof.* To prove correctness, it suffices to show that $k(h) = k(g)$ and that there does not exist $h \in h$ such that $k(h \setminus \{h\}) = k(h)$. Initially, we have $k(h) = k(g)$, and the field generated by $h$ does not change during the iterations. Thus, it suffices to show that such $h$ does not exist. Assume that it does. Then, $h \in k(h \setminus \{h\})$. However, at some iteration $i \in \{1, \dots, m\}$, we had $h = h_i$, so $h$ would have been removed from $h$, which is a contradiction.

There are $m$ independent applications of Algorithm 4, so we combine Proposition 2.21 with the Bernoulli inequality to obtain that the overall probability of correctness is at least

$$ \left(1 - \frac{\varepsilon}{m}\right)^m \ge 1 - \varepsilon. $$

□

**Remark 3.2 (On simplicity).** We note that result depends on the order of elements in the input $g$. We may order $g$ so that simple elements come last, so that they are used to remove other elements from the list. In practice, we use the following total order. Let $a/b, p/q \in k(x)$. Recall that for $f \in k[x]$, $d_f$ is the degree of $f$ and $s_f$ is the number of terms in $f$. We will ensure $d_a \ge d_b$ and $d_p \ge d_q$ by taking the reciprocal of the functions. We say that $a/b$ is simpler than $p/q$:

* if $d_a + d_b < d_p + d_q$, or, in the case of equality,

* if $s_a + s_b < s_p + s_q$, or, in the case of equality,

* if $d_b < d_q$, or, in the case of equality,

* if the largest monomial present in only one of $a$ or $p$ (when such exists) appears in $p$.
## 3.2 Low-degree coefficients of Gröbner bases over rational function fields

As we have seen, simple generators can be obtained from the low-degree coefficients of Gröbner bases of OMS ideals. In this section, we show how to compute those coefficients efficiently using sparse interpolation. We note that our algorithm is not restricted to OMS ideals and can be used with arbitrary ideals over the field of rational functions.

Let $f_1, \dots, f_m \in k(\mathbf{x})[\mathbf{y}]$, where $\mathbf{x} = (x_1, \dots, x_n)$ and $\mathbf{y} = (y_1, \dots, y_l)$. In this section, all computations are done with respect to the same fixed monomial ordering. Let $G$ be the reduced Gröbner basis of $\langle f_1, \dots, f_m \rangle$. We have $G = (g_1, \dots, g_\eta) \subset k(\mathbf{x})[\mathbf{y}]$, where for every $i = 1, \dots, \eta$,

$$g_i = \sum_{j=1}^{s_{g_i}} c_{i,j} \Pi(\mathbf{y}^{\alpha_{i,j}}),$$

where for every $j = 1, \dots, s_{g_i}$ we have $c_{i,j} = p_{i,j}/q_{i,j} \in k(\mathbf{x}) \setminus \{0\}$ with coprime $p_{i,j}, q_{i,j} \in k[\mathbf{x}] \setminus \{0\}$, and $\alpha_{i,j} \in \mathbb{N}_0^l$ are pairwise different.

Let $d$ be an integer and

$$ \begin{aligned} A_{\le d} &:= \{(i,j) \mid \deg p_{i,j} + \deg q_{i,j} \le d, i = 1, \dots, \eta, j = 1, \dots, s_{g_i}\}, \\ C_{\le d} &:= \{c_{i,j} \mid (i,j) \in A_{\le d}\} \subset k(\mathbf{x}). \end{aligned} \quad (5) $$

Note $C_{\le d}$ is a subset of the coefficients of $G$. We wish to compute $C_{\le d}$ without computing the whole $G$.

**Lemma 3.3.** For every $f_1, \dots, f_m \in k(\mathbf{x})[\mathbf{y}]$ and any monomial ordering, there exists a nonzero polynomial $g(\mathbf{x}) \in k[\mathbf{x}]$ with the following property. For every $\mathbf{a} \in k^n$ such that $g(\mathbf{a}) \ne 0$ the sets of polynomials in $k[\mathbf{y}]$ are well-defined and equal:

* the reduced Gröbner basis of $\langle f_1(\mathbf{a})(\mathbf{y}), \dots, f_m(\mathbf{a})(\mathbf{y}) \rangle$;

* the image of the reduced Gröbner basis of $\langle f_1, \dots, f_m \rangle \subset k(\mathbf{x})[\mathbf{y}]$ under the substitution $\mathbf{x} \to \mathbf{a}$.

*Proof.* Consider the run of the Buchberger's algorithm for $f_1, \dots, f_m$ over $k(\mathbf{x})$. We set $g$ to be the least common multiple of all the denominators occurring in the input and during the computation as well as the numerators of the leading coefficients of all input and intermediate polynomials. Consider $\mathbf{a} \in k^n$ such that $g(\mathbf{a}) \ne 0$. The Buchberger's algorithm, in addition to polynomial operations, queries the coefficients of polynomials in the following cases:

1. it queries the leading terms to form critical pairs and compute S-polynomials;

2. it queries all the nonzero coefficients of a polynomial when reducing it.

We consider the runs of the algorithm for $f_1(\mathbf{y}), \dots, f_m(\mathbf{y})$ and for its specialization $f_1(\mathbf{a})(\mathbf{y}), \dots, f_m(\mathbf{a})(\mathbf{y})$. Due to the assumption $g(\mathbf{a}) \ne 0$, the queries of the first type will yield the same leading monomials and the returned coefficients will commute with the specialization. For the queries of the second type, the only possible difference is if a non-leading coefficient has vanished under specialization. In this case, we will assume that its coefficient is still returned but is equal to zero, so there will be an extra reduction performed with zero coefficient which will not affect the result. Therefore, the polynomial operations performed will be the same regardless if the specialization was performed before the Gröbner basis computation or after. $\square$

Let $\text{coeff}(G, i, \alpha)$ denote the coefficient in front of the monomial $\Pi(\mathbf{y}^\alpha)$ in the $i$-th polynomial in $G$, or zero if the $i$-th polynomial or such monomial are not present in $G$; we shall reuse this notation for different $G$.

Let $i \in \{1, \dots, \eta\}$ and $j \in \{1, \dots, s_{g_i}\}$. Let $\text{BB}_{i,j}$ be a blackbox that, given as input an arbitrary $\mathbf{a} \in k^n$, computes $G_{\mathbf{a}} \subset k[\mathbf{y}]$, the reduced Gröbner basis of $\langle f_1(\mathbf{a})(\mathbf{y}), \dots, f_m(\mathbf{a})(\mathbf{y}) \rangle$, and returns $\text{coeff}(G_{\mathbf{a}}, i, \alpha_{i,j})$. More precisely, three outcomes are possible:

* $\text{BB}_{i,j}(\mathbf{a})$ returns `FAIL`. This happens when a denominator in $f$ vanishes under $\mathbf{x} \to \mathbf{a}$.

* $\text{BB}_{i,j}(\mathbf{a}) \ne c_{i,j}(\mathbf{a})$. This may happen when the image of $G$ under $\mathbf{x} \to \mathbf{a}$ and $G_{\mathbf{a}}$ do not coincide.

* $\text{BB}_{i,j}(\mathbf{a}) = c_{i,j}(\mathbf{a})$. This happens when none of the above happens. By Lemma 3.3, there exists a nonzero polynomial $g \in k[\mathbf{x}]$, such that for every $\mathbf{a} \in k^n$ if $g(\mathbf{a}) \ne 0$ then $\text{BB}_{i,j}(\mathbf{a}) = c_{i,j}(\mathbf{a})$.

The idea of the algorithm is as follows. First, consider the case when we have access to the support of $G$, that is, $\eta$ is known and for every $i = 1, \dots, \eta$ and $j = 1, \dots, s_{g_i}$ we have access to $\alpha_{i,j}$. Let $S \subseteq k \setminus \{0\}$. The output of the blackboxes is normalized consistently across the evaluations, since the Gröbner bases are reduced. So, for every $i = 1, \dots, \eta$ and $j = 1, \dots, s_i$, we can compute $\text{deg} p_{i,j}$ and $\text{deg} q_{i,j}$ by
applying Algorithm 3 to $BB_{i,j}$, $d$, and $S$. Consequently, this yields $A_{\le d}$. Let $\omega \in (k \setminus \{0\})^{n+1}$ be an admissible ratio as in Section 2.4. Then, for every $(i, j) \in A_{\le d}$, we may pass $(BB_{i,j}, S, \deg p_{i,j}, \deg q_{i,j}, \omega)$ to sparse interpolation from Algorithm 2 to recover the sparse representation of $c_{i,j} \in C_{\le d}$.

Now, for every $i = 1, \dots, \eta$ and $j = 1, \dots, s_{g_i}$, we can obtain $\alpha_{i,j}$ with high probability by computing $G_\sigma$ for $\sigma \in k^n$ whose components are chosen independently and uniformly from $S$. This gives rise to Algorithm 6.

<caption>**Algorithm 6 Computing low-degree coefficients of Gröbner bases**</caption>

**Input:** $f_1, \dots, f_m \in k(\mathbf{x})[\mathbf{y}]$; an integer $d$; a subset $S \subseteq k \setminus \{0\}$.

**Output:** $C_{\le d} \subset k(\mathbf{x})$ as in (5), or FAIL.

(Step 1) Pick an admissible ratio $\omega \in (k \setminus \{0\})^{n+1}$, as in Section 2.4; choose the coordinates of $\sigma \in k^n$ independently and uniformly from $S \subseteq k \setminus \{0\}$.

(Step 2) For every $i = 1, \dots, \eta$ and every $j = 1, \dots, s_{g_i}$, obtain $\alpha_{i,j} \in \mathbb{N}_0^t$ by computing $G_\sigma$, the reduced Gröbner basis of $\langle f_1(\sigma)(\mathbf{y}), \dots, f_m(\sigma)(\mathbf{y}) \rangle \subseteq k[\mathbf{y}]$.

(Step 3) For every $i = 1, \dots, \eta$ and every $j = 1, \dots, s_{g_i}$, construct the blackbox $BB_{i,j}$ that accepts as input $a \in k^n$, computes $G_a$, the reduced Gröbner basis of $\langle f_1(\mathbf{a})(\mathbf{y}), \dots, f_m(\mathbf{a})(\mathbf{y}) \rangle \subseteq k[\mathbf{y}]$, and returns either $\text{coeff}(G_a, i, \alpha_{i,j})$ or FAIL.

(Step 4) For every $i = 1, \dots, \eta$ and every $j = 1, \dots, s_{g_i}$, compute $\deg p_{i,j}$ and $\deg q_{i,j}$, or STOPPED if $\deg p_{i,j} + \deg q_{i,j} > d$, by applying Algorithm 3 to $(BB_{i,j}, d, S)$. If any of the applications failed, return FAIL.

(Step 5) Deduce $A_{\le d}$ as in (5).

(Step 6) For every $(i, j) \in A_{\le d}$, apply Algorithm 2 to $(BB_{i,j}, S, \deg p_{i,j}, \deg q_{i,j}, \omega)$ to obtain $p_{i,j}/q_{i,j} \in C_{\le d}$. If any of the applications failed, return FAIL.

(Step 7) Return $C_{\le d}$.

**Proposition 3.4.** Let $h := \deg g$ as in Lemma 3.3. Let $\nu = s_{g_1} + \dots + s_{g_\eta}$. Let $D$ be the maximum of $\deg p_{i,j} + \deg q_{i,j} + 2$ over all $i = 1, \dots, \eta$ and $j = 1, \dots, s_{g_i}$. Let $s$ be the maximum of $\max(s_{p_{i,j}}, s_{q_{i,j}})$ over all $(i, j) \in A_{\le d}$. The probability that Algorithm 6 succeeds and produces a correct result is at least

$$ \underbrace{\left(1 - \frac{h}{|S|}\right)}_{\text{(Step 2)}} \underbrace{\left(1 - \frac{\nu D(2 + 3h + 3D)}{|S|}\right)}_{\text{(Step 4)}} \underbrace{\left(1 - \frac{\nu D(\lceil \log_2 s \rceil + 8s + 4sh + 8sD)}{|S|}\right)}_{\text{(Step 6)}} $$

*Proof.* The blackboxes constructed in (Step 3) might return a wrong result and thus, technically, do not conform to the specifications of Algorithms 2 and 3. Our goal is to reconcile them. First we observe that Algorithm 6 always terminates because Algorithms 2 and 3 always terminate.

Next, for the sake of theoretical argument, we will consider a modified version of Algorithm 6 which we will call Algorithm 6*. The difference between these algorithms is that the blackbox constructed in (Step 3) of Algorithm 6* will check if the computed Gröbner basis is equal to the specialization of the Gröbner basis over $k(\mathbf{x})$ and return FAIL if it is not. This check can be performed, e.g., by explicitly computing the Gröbner basis over $k(\mathbf{x})$ (recall that Algorithm 6* is a purely theoretical construction). This makes the constructed blackboxes comply with the requirements of Algorithms 2 and 3, and the set of inputs yielding failure is still contained in the zero set of $g$ from Lemma 3.3. We note that, if Algorithm 6* returns a result, the same result would be returned by Algorithm 6 with the same input and the same random choices. Therefore, the probability that Algorithm 6 will return a correct result is at least the correctness probability for Algorithm 6*. Algorithm 6* succeeds and returns a correct result when (Step 2) computes the correct exponents and all of the applications of Algorithms 2 and 3 succeed and return correct results. Note:

* The computation in (Step 2) computes the correct exponents with probability at least $1 - h/|S|$.

* The individual applications of Algorithms 2 and 3 are pairwise independent, and all of those are independent from the computation in (Step 2).

So we apply Propositions 2.8 and 2.9 and the claim follows. $\square$

**Remark 3.5.** In Algorithm 6, it suffices to compute a single Gröbner basis to simultaneously obtain the evaluations of all blackboxes, since the admissible ratio is shared among the blackboxes. In particular, in the notation of Proposition 3.4, Algorithm 6 computes $\Theta sd$ Gröbner bases in $k[\mathbf{y}]$ to interpolate the result.
### 3.3 Finding polynomial generators

One specific way to find simple generators is to search for polynomials of low degree in the field $\mathcal{E}$. Polynomials in a subfield of $k(\mathbf{x})$ can be found using the following simple lemma.

**Lemma 3.6.** Let $g_1, \dots, g_m \in k(\mathbf{x})$ and $p \in k[\mathbf{x}]$. Then $p \in k(\mathbf{g})$ is and only if $p(\mathbf{y})$ is equivalent to an element of $k(\mathbf{x})$ modulo eOMS$_g$.

*Proof.* If $p \in k(\mathbf{g})$, then $p(\mathbf{y}) - p(\mathbf{x}) \in \text{eOMS}_g$ by Lemma 2.12. In the opposite direction, assume that there is $q(\mathbf{x}) \in k(\mathbf{x})$ such that $p(\mathbf{y}) - q(\mathbf{x}) \in \text{eOMS}_g$. By the definition of eOMS$_g$, every element of eOMS$_g$ vanishes under substitution $\mathbf{y} \to \mathbf{x}, t \to Q(\mathbf{x})$. Then $q(\mathbf{x}) = p(\mathbf{x})$, so $p \in k(\mathbf{g})$ by Lemma 2.12. $\square$

**Notation 3.7.** Let $I \subset k[\mathbf{x}]$ be an ideal and $p \in k[\mathbf{x}]$ be a polynomial. We fix the degrevlex monomial ordering on $k[\mathbf{x}]$; another ordering could be used as well. By $\text{NF}_I^+(p)$ we will denote the image of the normal form of $p$ with respect to the Gröbner basis of $I$ in the quotient space $k[\mathbf{x}]/k \cdot 1$ (where we factor a $k$-vector space $k[\mathbf{x}]$ by a subspace spanned by $1 \in k[\mathbf{x}]$). In other words, $\text{NF}_I^+(p)$ is the normal form of $p$ with removed constant term.

Using this notation, Lemma 3.6 can be restated as $p \in k(\mathbf{g}) \iff \text{NF}_{\text{eOMS}_g}^+(p) = 0$.

**Algorithm 7 Finding low-degree polynomial generators (cf. [69, Section 3.3])**

**Input:** $g_1, \dots, g_m \in k(\mathbf{x})$, with char $k = 0$; integer $\delta$ and subset $S \subseteq k \setminus \{0\}$;

**Output:** a basis of the vector space $k(\mathbf{g}) \cap \{p(\mathbf{x}) \in k[\mathbf{x}] | \text{deg } p \le \delta\}$.

(Step 1) Set $V$ to be the $k$-vector space spanned by the monomials in $\mathbf{y} = (y_1, \dots, y_n)$ of degree at most $\delta$.

(Step 2) Set $s := \text{dim } V$.

(Step 3) Repeat

(a) Set $\mathbf{a}$ be a vector in $k^n$ with the coordinates sampled uniformly at random from $S$.

(b) Set $V$ to be the kernel of the linear map $\text{NF}_{\text{eOMS}_g(a)}^+(V) : V \to k[\mathbf{y}, t]/k \cdot 1$.

The range of the map belongs to the (finite-dimensional) subspace of polynomials of $\text{deg} \le \delta$

(c) If $\text{dim } V = s$, then **return** a basis of $V$ substituting $\mathbf{y} \to \mathbf{x}$.

(d) Set $s := \text{dim } V$.

**Proposition 3.8.** Algorithm 7 terminates. For every $g_1, \dots, g_m, \delta$, and $\varepsilon > 0$, there exists $D_0$ such that, if $|S| > D_0$, then the result is correct with probability at least $1 - \varepsilon$.

**Lemma 3.9.** Let $M(\mathbf{x})$ be a matrix over $k(\mathbf{x})$, and let $r$ be the dimension of the $k$-vector space spanned by the rows of $M$. Consider a matrix $M(\mathbf{x}_1, \dots, \mathbf{x}_r) := (M(\mathbf{x}_1) \dots M(\mathbf{x}_r))$ obtained by concatenating matrices $M(\mathbf{x}_1), \dots, M(\mathbf{x}_r)$, where $\mathbf{x}_1, \dots, \mathbf{x}_r$ are independent tuples of indeterminates. Then $M(\mathbf{x}_1, \dots, \mathbf{x}_r)$ has rank $r$.

Furthermore, for every $s < r$

$$ \text{rank } M(\mathbf{x}_1, \dots, \mathbf{x}_s) = \text{rank } M(\mathbf{x}_1, \dots, \mathbf{x}_{s+1}) \implies \text{rank } M(\mathbf{x}_1, \dots, \mathbf{x}_s) = r. $$

*Proof.* Since every $k$-linear dependence between the rows of $M(\mathbf{x})$ yields a $k$-linear dependence between the rows of $M(\mathbf{x}_1, \dots, \mathbf{x}_r)$, the rank of the latter is at most $r$. Let $c$ be the number of columns in $M(\mathbf{x})$. We introduce a new variable $y$ and consider a vector $v(\mathbf{x}, y) := M(\mathbf{x}) \cdot (1 \ y \ \dots \ y^{c-1})^T$. The $k$-linear span of the entries of the vector has also dimension $r$. Therefore, by [98, Lemma 19], the matrix $(v(\mathbf{x}_1, y_1) \ \dots \ v(\mathbf{x}_r, y_r))$ has rank $r$. Since this matrix is a product of $M(\mathbf{x}_1, \dots, \mathbf{x}_r)$ and a block matrix with blocks of the form $(1 \ y_i \ \dots \ y_i^{c-1})^T$, the rank of $M(\mathbf{x}_1, \dots, \mathbf{x}_r)$ is at least $r$. This proves rank $M(\mathbf{x}_1, \dots, \mathbf{x}_r) = r$.

To prove the second claim, we observe that the equality rank $M(\mathbf{x}_1, \dots, \mathbf{x}_s) = \text{rank } M(\mathbf{x}_1, \dots, \mathbf{x}_{s+1})$ implies that the columns of $M(\mathbf{x}_{s+1})$ belong to the column space of $M(\mathbf{x}_1, \dots, \mathbf{x}_s)$. Then the same will be true for $M(\mathbf{x}_{s+2}), \dots, M(\mathbf{x}_r)$, so rank $M(\mathbf{x}_1, \dots, \mathbf{x}_s) = \text{rank } M(\mathbf{x}_1, \dots, \mathbf{x}_r) = r$. $\square$

**Proof of Proposition 3.8.** For proving termination, we observe that, at each iteration of the loop on (Step 3), either the dimension of $V$ becomes smaller or the algorithm returns. Thus, the algorithm will return after a finite number of iterations of the loop.
We apply Lemma 3.3 to the generators of the ideal eOMS<sub>g</sub> and denote the resulting polynomial by $b(\mathbf{x}) \in k[\mathbf{x}]$. Let $m_1, \dots, m_N$ be the monomials of degree at most $\delta$ in $\mathbf{y}$. Then $V$ at (Step 1) is spanned by $m_1, \dots, m_N$. Let $N_1$ be the dimension of the $k(\mathbf{x})$-space of polynomials of degree at most $\delta$ in $\mathbf{y}$ and $t$ without constant term. We choose a $k(\mathbf{x})$-basis in this space, then we can write $\text{NF}_{\text{eOMS}_g}^{+}(m_1), \dots, \text{NF}_{\text{eOMS}_g}^{+}(m_N)$ in the basis as vectors in $k(\mathbf{x})^{N_1}$. We consider an $N \times N_1$-matrix $M(\mathbf{x})$ whose rows are these vectors. By Lemma 3.6, polynomials of degree at most $\delta$ in $k(\mathbf{g})$ are in one-to-one correspondence with the $k$-linear dependencies between the rows of $M(\mathbf{x})$. We denote by $r$ the dimension of the $k$-vector space of these dependencies. By Lemma 3.9, matrix $M(\mathbf{x}_1, \dots, \mathbf{x}_r)$ has rank $r$. Let $\ell$ be the smallest integer such that $M(\mathbf{x}_1, \dots, \mathbf{x}_\ell)$ has rank $r$. Let the ranks of $M(\mathbf{x}_1), M(\mathbf{x}_1, \mathbf{x}_2), \dots, M(\mathbf{x}_1, \dots, \mathbf{x}_\ell)$ be $r_1, r_2, \dots, r_\ell = r$, respectively. By the second part of Lemma 3.9, $r_1 < r_2 < \dots < r_\ell$. For every $1 \le i \le \ell$, let $q_i(\mathbf{x}_1, \dots, \mathbf{x}_i)$ be the determinant of any nonsingular $r_i \times r_i$ minor of $M(\mathbf{x}_1, \dots, \mathbf{x}_i)$. We denote by $q(\mathbf{x}_1, \dots, \mathbf{x}_\ell)$ the product $q_1 \cdots q_\ell$.

Let $\mathbf{a}_1, \mathbf{a}_2, \dots$ be the integer vectors to be generated in the subsequent iterations of the loop from (Step 3) (we consider an infinite sequence regardless of the actual number of iterations). We will prove that if

$$b(\mathbf{a}_1) \dots b(\mathbf{a}_\ell) b(\mathbf{a}_{\ell+1}) q(\mathbf{a}_1, \dots, \mathbf{a}_\ell) \neq 0, \qquad (6)$$

then the algorithm will return correct result at the $(\ell+1)$-st iteration of the loop at (Step 3). Consider the $i$-th iteration of the loop. We will denote the values of $V$ and $s$ at the beginning of the iteration by $V_i$ and $s_i$, respectively. Since $b(\mathbf{a}_i) \neq 0$, computation of the normal form with respect to eOMS<sub>g</sub> and evaluation at $\mathbf{a}_i$ commute. Hence, the kernel of the map $\text{NF}_{\text{eOMS}_g}^{+}(\mathbf{a}_i)$ restricted to the polynomials of degree $\le \delta$ will be in bijective correspondence with $k$-linear relations between the rows of $M(\mathbf{a}_i)$. Therefore, $V_{i+1}$ will be isomorphic to the space of all $k$-linear relations between the rows of $M(\mathbf{a}_1, \dots, \mathbf{a}_i)$. Since $q_i(\mathbf{a}_1, \dots, \mathbf{a}_i) \neq 0$, we will have $s_{i+1} = r_i$. Thus, the algorithm will not terminate at the first $\ell$ iterations but will terminate at the $(\ell+1)$-st. Furthermore, since the elements of the resulting $V$ will be in correspondence with the $k$-relations on the rows of $M(\mathbf{a}_1, \dots, \mathbf{a}_{\ell+1})$, the final $V$ will contain the space of polynomials of degree $\le \delta$ in $k(\mathbf{g})$. Since $\dim V = r$, this containment will be in fact an equality.

Thus, if the condition (6) is fulfilled, the algorithm will return a correct result. By the Schwarz-Zippel lemma (Lemma 2.1), the probability of this can be made arbitrary close to 1 by increasing the size of $S$. $\square$

**Remark 3.10 (Modular version (cf. Remark 2.23)).** In practice, the computation will be performed modulo a prime number $p$ (and the coordinates for $\mathbf{a}$'s will be sampled from the whole $\mathbb{F}_p$). The proof of Proposition 3.8 implies that, for any given $\varepsilon > 0$, the probability of correctness will be at least $1 - \varepsilon$ for large enough $p$.

**Remark 3.11 (On computation of $k(\mathbf{g}) \cap k[\mathbf{x}]$).** Polynomials in $k(\mathbf{g})$ form not only a linear subspace but also a subalgebra. It is therefore tempting not to restrict ourselves to fixed $\delta$ but compute generators of this subalgebra. However, as the negative solution to the Hilbert 14-th problem [102] shows, this subalgebra may be not finitely generated (even in the three-variable case [83]).

# 4 Main algorithm

## 4.1 Algorithm

In this section, we put the subroutines described above together to present our main algorithm.
**Algorithm 8 Finding simple generators of a rational function field**

**Input:** $\mathbf{g} = (g_1, \dots, g_m) \subset k(\mathbf{x})$, where $\text{char } k = 0$ and we denote $\mathcal{E} := k(\mathbf{g})$; a real number $0 < \varepsilon < 1$.

**Output:** $\mathbf{h} = (h_1, \dots, h_\ell) \subset k(\mathbf{x})$, such that $k(\mathbf{h}) = k(\mathbf{g})$ (and $\mathbf{h}$ is potentially simpler than $\mathbf{g}$).

(Step 1) Set $d := 1$ and $\tau := 1$ and choose a finite subset $S \subset k \setminus \{0\}$.

(Step 2) Compute $c_1, \dots, c_M \in k(\mathbf{x})$, the coefficients up to degree $d$ of the Gröbner basis of $\text{eOMS}_g$, by applying Algorithm 6 to the generators of $\text{eOMS}_g$ and $d, S$.

(Step 3) Check if $\mathcal{E} = k(c_1, \dots, c_M)$ with probability at least $1 - \frac{\varepsilon}{\tau}$ using Algorithm 4 and Remark 2.25. If the equality does not hold, set $d := 2d$ and **go to (Step 2)**.

(Step 4) Compute $m_1, \dots, m_N \in k(\mathbf{x})$, polynomials up to degree $\delta$ that belong to $\mathcal{E}$, by applying Algorithm 7 to $c_1, \dots, c_M$ with $\delta$ set to a small constant, say 2 or 3, and the set $S$.

(Step 5) Sort the functions $(g_1, \dots, g_m, c_1, \dots, c_M, m_1, \dots, m_N)$ using a heuristic simplicity measure and obtain a list $(h_1, \dots, h_{m+N+M})$; for example, as in Remark 3.2.

(Step 6) For each $i = 1, \dots, m+N+M$, check if $h_i \in k(h_1, \dots, h_{i-1})$ with probability at least $1 - \frac{\varepsilon}{\tau}$ using Algorithm 4. Denote the list of the functions for which the inclusion did not hold by $(\hat{h}_1, \dots, \hat{h}_s)$. Optionally, also minimize by applying Algorithm 5 to $(\hat{h}_s, \dots, \hat{h}_1)$.

(Step 7) Check if $\mathcal{E} = k(\hat{h}_1, \dots, \hat{h}_s)$ with probability at least $1 - \frac{\varepsilon}{2\tau}$ using Algorithm 4 and Remark 2.25.

(Step 8) If the check has passed, **return h** := $(\hat{h}_1, \dots, \hat{h}_s)$. Otherwise, choose a larger subset $S \subset k \setminus \{0\}$, set $\tau := \tau + 1$, and **go to (Step 2)**.

**Proposition 4.1.** *Algorithm 8 terminates with probability one. The output of Algorithm 8 is correct with probability at least $1 - \varepsilon$.*

*Proof.* We prove the probability of termination. The algorithm does not terminate if and only if the check in (Step 7) never passes. The probability of that is $p_1 \cdot p_2 \cdot \dots$, where at the $\tau$-th iteration of the algorithm $p_\tau$ is the probability that (Step 7) does not pass. To show that this product is zero, it suffices to show that there exists a positive integer $n$ such that for every $\tau > n$ we have $p_\tau < 0.42$. We denote by $e_0$ the smallest positive integer such that the coefficients up to degree $2^{e_0}$ of the Gröbner basis of $\text{eOMS}_g$ generate the whole $\mathcal{E}$. Then, we fix these coefficients (and denote their number by $M_0$), and let $\rho(|S|)$ be the probability that Algorithm 7 applied to these coefficients returns a correct result. Let $N_0$ be the number of polynomial generators in this correct result. We fix a positive integer $\tau$ and consider the $\tau$-th iteration of the algorithm (in other words, the pass from (Step 2) to (Step 7) preceding the $\tau$-th check in (Step 7).) Consider the following events:

* The applications of Algorithm 6 and Algorithm 4 in (Step 2) and (Step 3) return correct result. In this case, there will be $e_0 + 1$ applications of these algorithms. The probability of correctness of applications of Algorithm 6 is of the form $f(|S|)$ where $f(|S|) \to 1$ when $|S| \to \infty$. Since $S$ becomes larger at each step, this probability tends to 1 as $\tau \to \infty$. The probability of correctness of the applications of Algorithm 4 is at least $1 - \frac{e_0\varepsilon}{\tau}$, which also tends to 1 as $\tau \to \infty$.

* The application of Algorithm 7 at (Step 4) is correct. This happens with probability $\rho(|S|) \to 1$ as $\tau \to \infty$ by Proposition 3.8.

* All the applications of Algorithm 4 in (Step 6) return correct result. This happens with probability at least $1 - \frac{(m+M_0+N_0)\varepsilon}{\tau}$. This value tends to 1 as $\tau \to \infty$.

* The check at (Step 7) returns correct result.

If all these events occur, then the check at (Step 7) passes. By the probability calculations above, the probability of this tends to 1 as $\tau$ approaches infinity. This means that $p_\tau$ tends to 0 as $\tau \to \infty$. In particular, there exists $n$ such that $p_\tau < 0.42$ for any $\tau > n$.

We prove the probability of correctness. At the $\tau$-th iteration of the algorithm, let $B_\tau$ denote the event that $\mathcal{E} = k(\hat{h}_1, \dots, \hat{h}_s)$, and let $E_\tau$ denote the event that the check in (Step 7) passes (that is, it reports that
$\mathcal{E} = k(\hat{h}_1, \dots, \hat{h}_s))$. At the $\tau$-th iteration, we have:

$$ 
\begin{align*}
P(\text{(Step 7) is correct}) &= P(E_{\tau} | B_{\tau})P(B_{\tau}) + P(\overline{E_{\tau}} | \overline{B_{\tau}})P(\overline{B_{\tau}}) &\ge 1 - \frac{\varepsilon}{2^{\tau}} \\
&\Rightarrow (1 - P(\overline{E_{\tau}} | B_{\tau}))P(B_{\tau}) + (1 - P(E_{\tau} | \overline{B_{\tau}}))P(\overline{B_{\tau}}) &\ge 1 - \frac{\varepsilon}{2^{\tau}} \\
&\Rightarrow P(\overline{E_{\tau}} | B_{\tau})P(B_{\tau}) + P(E_{\tau} | \overline{B_{\tau}})P(\overline{B_{\tau}}) &\le \frac{\varepsilon}{2^{\tau}} \\
&\Rightarrow P(E_{\tau} | \overline{B_{\tau}})P(\overline{B_{\tau}}) &\le \frac{\varepsilon}{2^{\tau}} \\
&\Rightarrow P(E_{\tau} \cap \overline{B_{\tau}}) &\le \frac{\varepsilon}{2^{\tau}}
\end{align*}
$$

Note that, at the $\tau$-th iteration, $E_{\tau} \cap \overline{B_{\tau}}$ is the event that an incorrect result is returned from the algorithm.

Let $T_{\tau}$ denote the event that the algorithm terminates at the $\tau$-th iteration with an incorrect result and that it did not terminate at every previous iteration. Then,

$$ P(T_{\tau}) \le P(E_{\tau} \cap \overline{B_{\tau}}) \le \frac{\varepsilon}{2^{\tau}}. $$

Since iterations of the algorithm are independent, the probability of an incorrect result is

$$ T_1 + T_2 + \dots \le \frac{\varepsilon}{2} + \frac{\varepsilon}{4} + \dots = \varepsilon \left( \frac{1}{2} + \frac{1}{4} + \dots \right) \le \varepsilon. $$

Since the algorithm terminates with probability one, the probability of a correct result and the probability of an incorrect result sum up to one. Therefore, the probability of a correct result is at least $1 - \varepsilon$. $\square$

**Remark 4.2** (Retaining the original generators). In (Step 5), we extend the list of candidate functions with the original $g_1, \dots, g_m$. Roughly speaking, this ensures that the output of Algorithm 8 is not worse than the input in terms of being simple. For example, consider $g = (x^i + y^i + z^i + u^i + v^i, i = 1, \dots, 5)$. The output of our algorithm for such input (using the heuristic from Remark 3.2) is

$$ 
\begin{align*}
x + y + z + u + v, x^2 + y^2 + z^2 + u^2 + v^2, x^3 + y^3 + z^3 + u^3 + v^3, \\
xyzu + xyzv + xyuv + xzuv + yzuv, xyzuv.
\end{align*}
$$

However, if we did not add the original generators to the list of candidates, the output would be less simple (both visually and according to the chosen heuristic):

$$ 
\begin{align*}
x + y + z + u + v, xy + xz + xu + xv + yz + yu + yv + zu + zv + uv, \\
xyz + xyu + xyv + xzu + xzv + xuv + yzu + yzv + yuv + zuv, \\
xyzu + xyzv + xyuv + xzuv + yzuv, xyzuv.
\end{align*}
$$

**Remark 4.3** (Modular version for $k = \mathbb{Q}$). In practice, we run the steps (Step 1)–(Step 6) of Algorithm 8 modulo a large prime number. Only the final check in (Step 7) is run over $k = \mathbb{Q}$, and its purpose is to guarantee the probability of correctness, as Proposition 4.1 would still hold in this case. Thus just prior to the final check we reconstruct the intermediate results using rational number reconstruction [126, Section 5.10], [127]. The use of Algorithms 4, 5 and 7 modulo a prime is justified due to Remarks 2.23 and 3.10.

**Remark 4.4** (Lazy specialization of eOMS$_g$). In practice, we do not explicitly construct the generators of eOMS$_g$, since they may be large. Indeed, let $p/q \in g$ with coprime $p, q \in k[x]$. The corresponding generator of eOMS$_g$ is $p(y)q(x) - q(y)p(x)$ and it has $O(s_p s_q)$ terms when considered as a polynomial in $k[x, y]$. In comparison, the corresponding generator of eOMS$_g(a)$ for $a \in k^n$ has only $O(s_p + s_q)$ terms. Therefore, whenever we need to compute a Gröbner basis of eOMS$_g(a)$ for $a \in k^n$, we construct the generators of eOMS$_g(a)$ on the fly by plugging $x \to a$ in (4).

**Remark 4.5** (Computing with zero-dimensional ideals). Field membership (Algorithm 4) specializes some indeterminates in the OMS ideal until the ideal becomes zero-dimensional. In theory, a similar specialization could be used for computing the coefficients of Gröbner bases in Algorithm 6. In this case, in some polynomials in the resulting ideal, the monomials may glue together, and their coefficients (which are multivariate rational functions) combined may become large. In our experiments, we observed that this makes interpolation harder, and overall slows down the algorithm despite using zero-dimensional ideals.

**Remark 4.6.** After we have computed the coefficients of the Gröbner basis that generate the field, we can use them instead of the original generators in the subsequent steps of the algorithm (as in (Step 4)). This typically speeds up the computation when these coefficients are simpler than the original generators.
**Remark 4.7** (Achieving minimality). The generating set computed by Algorithm 8 may be not minimal with respect to inclusion (see Example 7.3) despite the filtering performed at (Step 6). Since the filtering checks membership against the subfield generated by the elements with smaller indices only, there could still be a redundant element at a smaller index that can be expressed in terms of all other generators. The minimality with respect to inclusion can be achieved using Algorithm 5. However, in some cases, the larger set may be preferable, see the discussion in Example 7.3.

## 4.2 Example

In this section, we run Algorithm 8 on an example and trace through the steps of the algorithm. We consider an example from the area of structural identifiability of dynamical models (Section 7.1) with the code name SEIR34, collected in [90]. The task consists of finding a simple generating set of the subfield $\mathcal{E} \subseteq \mathbb{Q}(k, N, \beta, \epsilon, \gamma, \mu, r)$ generated by

$$ \begin{gathered} N, \epsilon + \gamma + 3\mu, -\epsilon - \gamma - 2\mu, \epsilon\mu + \gamma\mu + 2\mu^2, \epsilon\gamma\mu + \epsilon\mu^2 + \gamma\mu^2 + \mu^3, \\ \frac{\beta r}{kN}, \frac{-\beta\epsilon r}{N}, \frac{\beta\epsilon r + \beta\gamma r + 2\beta\mu r}{kN}, \frac{\beta\epsilon\gamma r + \beta\epsilon\mu r + \beta\gamma\mu r + \beta\mu^2 r}{kN}. \end{gathered} $$

These nine functions are the input to our Algorithm 8. When executing the steps of the algorithm, we will suppose that all steps return a correct result. When computing Gröbner bases, we will be using the degrevlex monomial ordering with $t > k > N > \beta > \epsilon > \gamma > \mu > r$ (where $t$ comes from the Rabinowitsch trick).

In (Step 2), the algorithm computes the coefficients of degree at most $d = 1$ of the Gröbner basis of the OMS ideal and gets

$$ (c_1, \dots, c_M) := (\mu, \epsilon + \gamma, N). \qquad (7) $$

These do not generate $\mathcal{E}$, thus the check in (Step 3) does not pass. Therefore, we run (Step 2) again, now with $d = 2$, and obtain the same result as (7) but additionally with two coefficients of higher degree: $\epsilon\gamma$ and $\frac{k}{\gamma}$. The check in (Step 3) still fails, so we run again with $d = 4$ and discover an additional coefficient: $\frac{\beta r}{\gamma}$. At this point,

$$ (c_1, \dots, c_M) := \left( \mu, \epsilon + \gamma, N, \epsilon\gamma, \frac{k}{\gamma}, \frac{\beta r}{\gamma} \right). $$

Finally, the check signals that $\mathcal{E} = \mathbb{Q}(c_1, \dots, c_M)$.

Then, (Step 4) computes some linearly independent polynomials up to degree $\delta = 2$ that belong to $\mathcal{E}$:

$$ (m_1, \dots, m_N) := (\mu, \epsilon + \gamma, N, k\epsilon, \epsilon\gamma). $$

Then, in (Step 5), the original generators and $(c_1, \dots, c_M)$ and $(m_1, \dots, m_N)$ are combined in one list, and the list is ordered using the heuristic from Remark 3.2, which ensures that simpler elements come first. In this example, we omit the original generators, as adding them does not change the final result. After removing duplicates, we have

$$ (h_1, \dots, h_{N'+M'}) := \left( \mu, N, \epsilon + \gamma, \epsilon\gamma, k\epsilon, \frac{k}{\gamma}, \frac{\beta r}{\gamma} \right). $$

Note that $\mathcal{E} = \mathbb{Q}(h_1, \dots, h_{N'+M'})$. The purpose of (Step 6) is to pick the simplest elements from this list that would form the final result. To this end, the algorithm makes a pass over the list and removes $\frac{k}{\gamma}$ since it can be expressed in terms of the previous elements in the list. We obtain the simplified generating set (minimization from Algorithm 5 has no effect):

$$ (\hat{h}_1, \dots, \hat{h}_s) := \left( \mu, N, \epsilon + \gamma, \epsilon\gamma, k\epsilon, \frac{\beta r}{\gamma} \right). $$

These functions are the output of our algorithm.

This result requires both the coefficients of the Gröbner basis and some of the polynomial generators. In particular, the function $k\epsilon$ is not present among the coefficients of the Gröbner basis. In turn, the subfield cannot be generated by $m_1, \dots, m_N$ alone, thus computing some coefficients of the Gröbner basis was necessary.

**Remark 4.8** (Algorithm parameters). We can increase the effort committed to simplification by computing Gröbner bases of OMS ideals in several monomial orderings, as in Example 2.17, and by increasing the parameter $\delta$ of Algorithm 7, and then combining the results. In this section's example, if we used $\delta = 3$ in (Step 4), we would additionally find that $\beta\epsilon r \in \mathcal{E}$, and hence obtain a generating set that contains only polynomials: $\mu, N, \epsilon + \gamma, \epsilon\gamma, k\epsilon, \beta\epsilon r$.
# 5 Implementation and performance

## 5.1 Implementation

We have implemented Algorithm 8 and the sub-algorithms over the field of rational numbers $k = \mathbb{Q}$ in Julia language [12]. Our implementation has a convenient user interface and is freely available in the package RationalFunctionFields.jl at

https://github.com/pogudingleb/RationalFunctionFields.jl <fnref n="1" />

One essential component is the algorithm for interpolating the coefficients of a Gröbner basis over the rational function fields, which is implemented in the sub-package ParamPunPam.jl. We use Nemo.jl [40] for polynomial arithmetic and Groebner.jl [34] for computing Gröbner bases over integers modulo a prime.

The experiments in this paper were run on a Linux machine with 128 GB of RAM and 32 × 13th Gen i9-13900+ CPU. We used Julia 1.12.4, Singular 4.4.1, and Maple 2025.

## 5.2 Modular computation

Input and output in our implementation of the main algorithm are rational functions over the rationals. Per Remark 4.3, for better performance, intermediate computations are carried out over finite fields. The final correctness check in (Step 7) computes Gröbner basis and normal forms over the rationals. Typically, the size of the rational numbers in the output of our algorithm is small in our applications (Section 7); this is also observed in practice (see, e.g., Table 4 and Appendix A).

**Remark 5.1** (Reconciling with reality). In our implementation, the final correctness check in (Step 7) computes Gröbner basis over the rationals using a Monte-Carlo probabilistic algorithm based on multi-modular computation [5], for which no practical bound on the probability of correctness is known. For this reason, our implementation may not guarantee the probability of correctness stated in Proposition 4.1.

## 5.3 Sparse interpolation over finite fields

For interpolating sparse multivariate rational functions over finite fields, we implemented the algorithm by van der Hoeven and Lecerf [124] discussed in Section 2.5. We also implemented the algorithm by Cuyt and Lee [29] because, in theory, it may require fewer blackbox evaluations; see [124, Remark 3.]. In most of our examples, the algorithms used the same number of evaluations, so we prefer van der Hoeven and Lecerf’s algorithm because it was easier to implement efficiently.

For interpolating sparse multivariate polynomials over a finite field $\mathbb{F}_q$, we use the Ben-Or and Tiwari algorithm [10] with the admissible ratio $\omega = (2, 3, 5, \dots, p_n)$ and recover exponents via integer factoring (see Section 2.4). Therefore, to interpolate polynomials of degree $d$, we must have $p_n^d < q$. To facilitate efficient arithmetic in $\mathbb{F}_q$, we choose $q$ to be a prime that fits in 64 bits. For example, with $d = 8$, the largest number of variables such that $p_n^d < 2^{64}$ is $n = 54$, which has been sufficient in our use cases. An alternative would have been Kronecker substitution ([126, Section 8.4.], [4]), which requires, roughly, $d^n < q$. Our choice between the two is motivated by the fact that we usually have $n > d$.

We double the number of evaluations until interpolation succeeds. We found this approach easier to implement in an efficient way than the early termination strategies (as, for example, in [73, 75]), since we could rely on the existing Páde approximation and root-finding functionality from Nemo.jl. Thus, our implementation overestimates the required number of blackbox evaluations by at most a factor of two.

## 5.4 Computing small degree coefficients of Gröbner bases

An appealing property of our algorithm is that the degree of the interpolated coefficients of the Gröbner basis of the OMS ideal is well-controlled. In Algorithm 8, (Step 2) and (Step 3) ensure that the coefficients of high degree are not interpolated when they are not necessary to generate the field (as, e.g., in Example 2.16). Along with being more efficient, this method also helps with simplification by filtering out generators of high degrees before even computing them.

<footnote n="1">¹The URL refers to the most up-to-date version of the code. The code as it had been at the time of this publication can be accessed at git revision f7d4943aa7432bd10c206b59f63909aca52c909d at: https://github.com/pogudingleb/RationalFunctionFields.jl/tree/f7d4943aa7432bd10c206b59f63909aca52c909d. In particular, the directory paper contains the sources of our experiments and examples.</footnote>
<table><caption>Table 1: The number of blackbox evaluations required to interpolate the coefficients of Gröbner bases of OMS ideals, using two methods: computing the full Gröbner basis with Algorithm 6, and using Algorithm 6 with the same *d* as in (Step 2) and (Step 3) of Algorithm 8. We write n/a if the computation is not feasible.</caption><thead><tr><th rowspan="2">Example</th><th colspan="2">Input</th><th colspan="2">Full Gröbner basis</th><th colspan="2">(Step 2) and (Step 3)</th></tr><tr><th># vars.</th><th>degrees</th><th>degrees</th><th># evals.</th><th>degrees</th><th># evals.</th></tr></thead><tbody><tr><td>MAPK-5 [89]</td><td>22</td><td>8, 6</td><td>12, 0</td><td>15,401</td><td>1, 0</td><td>16</td></tr><tr><td>Pharm [32]</td><td>7</td><td>26, 5</td><td>2, 0</td><td>20</td><td>2, 0</td><td>20</td></tr><tr><td>Fujita [46]</td><td>16</td><td>9, 5</td><td>4, 7</td><td>129</td><td>1, 1</td><td>42</td></tr><tr><td>Bilirubin [94]</td><td>7</td><td>4, 0</td><td>4, 2</td><td>278</td><td>4, 2</td><td>278</td></tr><tr><td>EAIHRD [45]</td><td>10</td><td>22, 6</td><td>10, 9</td><td>n/a</td><td>4, 4</td><td>342</td></tr><tr><td>SLIQR [30]</td><td>6</td><td>6, 0</td><td>7, 6</td><td>502</td><td>7, 6</td><td>502</td></tr><tr><td>LinComp2 [1]</td><td>7</td><td>5, 0</td><td>381, 373</td><td>n/a</td><td>4, 2</td><td>4118</td></tr></tbody></table>

We will illustrate this claim using several examples arising from the problem of structural identifiability (see Section 7.1). Each example is given as a list of rational functions – the generating set to be simplified. We will consider two methods. For each example, we apply Algorithm 6 to the generators of the OMS ideal of the input functions, but the first method computes the full Gröbner basis (by specifying some large degree bound *d*), and the second method uses the degree *d* from (Step 2) and (Step 3) of Algorithm 8.

In Table 1, for each example, we report: the number of variables; the largest degree of the numerator and denominator among the input rational functions and the output rational functions, for both methods; and the number of evaluations required to interpolate the output, for both methods.

From Table 1, we may conclude that

* The maximal degrees in the output typically are not large (even in the full Gröbner bases), and often are smaller than those of the input. This supports the expectations that OMS ideals are useful for simplification.

* In some cases, the field can be generated by elements with degrees smaller than that of the full Gröbner basis, and exploiting this leads to fewer evaluations required for interpolation. For example, this makes the computation feasible for the EAIHRD and LinComp2 examples, and reduces the running time from 16 minutes to 2 seconds for the MAPK-5 example.

**Remark 5.2** (Lowering the degrees even further). The strategy of doubling the guess for the degree *d* until the check in (Step 3) passes may be suboptimal. For instance, in the SLIQR example in Table 1, the maximal total degrees of the numerator and denominator of a coefficient in the full Gröbner basis are 7 and 6, respectively; the degrees obtained in (Step 2) and (Step 3) are the same. However, the check in (Step 3) would pass for degrees as low as (5, 4). Thus, with the current strategy, we interpolate coefficients of higher degrees than necessary. This may significantly increase the number of blackbox evaluations as the maximum number of terms $\binom{n+d}{n}$ grows rapidly with the degree. Since the check in (Step 3) is cheap, we could search for a sharper guess for the degree, for example, by increasing *d* in arithmetic progression with a small step. The current strategy is a compromise with the fact that we do not cache blackbox evaluations across the calls to Algorithm 6, thus we would like to minimize the number of iterations.

## 5.5 Speeding up computation of many Gröbner bases

The efficiency of Algorithm 6 depends on the ability to evaluate the blackboxes many times, where each evaluation involves computing a Gröbner basis over a finite field. As we have seen, one might require many evaluations to interpolate the result.

We use Gröbner tracing, as described by Traverso [122], to speed up the computation of Gröbner bases. The presentation in the original paper is general, however, for completeness, we briefly present the main ideas while adapting them to our case.
<table><caption>Table 2: The time to compute a single Gröbner basis over integers modulo a prime using the F4 algorithm with and without tracing. We denote microseconds by µs.</caption><thead><tr><th>Example</th><th>F4 algorithm</th><th>F4 with tracing (Apply)</th><th>Speedup</th></tr></thead><tbody><tr><td>MAPK-5 [89]</td><td>39 ms</td><td>8 ms</td><td>4.9</td></tr><tr><td>Pharm [32]</td><td>1.91 s</td><td>0.36 s</td><td>5.3</td></tr><tr><td>Fujita [46]</td><td>1.89 ms</td><td>0.54 ms</td><td>3.5</td></tr><tr><td>Bilirubin [94]</td><td>334 µs</td><td>129 µs</td><td>2.6</td></tr><tr><td>EAIHRD [45]</td><td>680 ms</td><td>112 ms</td><td>6.1</td></tr><tr><td>SLIQR [30]</td><td>116 µs</td><td>39 µs</td><td>3.1</td></tr><tr><td>LinComp2 [1]</td><td>68 ms</td><td>37 ms</td><td>1.8</td></tr></tbody></table>

We describe the method for OMS ideals but the algorithm can be applied to any ideal over a rational function field. Let $g = (g_1, \dots, g_m) \subset k(\mathbf{x})$ with $\mathbf{x} = (x_1, \dots, x_n)$ and eOMS$_g \subset k(\mathbf{x})[\mathbf{y}]$ be its OMS ideal. Recall that eOMS$_g(\mathbf{a})$ denotes the specialization of the OMS ideal under $\mathbf{x} \to \mathbf{a}$, where $\mathbf{a} \in k^n$. Fix a monomial ordering on $\mathbf{y}$. Let $\mathbf{a}_1, \dots, \mathbf{a}_\ell \in k^n$. For $i = 1, \dots, \ell$, we wish to compute the reduced Gröbner basis of eOMS$_g(\mathbf{a}_i)$. With the use of tracing, the computation is organized in two steps:

(Learn) Run a modified Buchberger's algorithm to compute the Gröbner basis of eOMS$_g(\mathbf{a}_1)$. The modified algorithm additionally records the information about which S-polynomials reduced to zero; this information could be stored as a collection of indices of the corresponding iterations of the algorithm.

(Apply) For $i = 2, \dots, \ell$, run a (differently) modified Buchberger's algorithm to compute the Gröbner basis of eOMS$_g(\mathbf{a}_i)$. This version runs as usual but does not reduce the S-polynomials that correspond to the ones that reduced to zero during the (Learn) stage, since those are expected to reduce to zero.

Correctness of this method can be analyzed using Lemma 3.3. Indeed, there exists a nonzero polynomial $g \in k(\mathbf{x})$ such that for every $\mathbf{b}_1, \mathbf{b}_2 \in k^n$ if $g(\mathbf{b}_1) \neq 0$ and $g(\mathbf{b}_2) \neq 0$ then the Buchberger's algorithm applied to the generators of eOMS$_g(\mathbf{b}_1)$ and eOMS$_g(\mathbf{b}_2)$ produces the same S-polynomials and the same reductions of S-polynomials (that is, with the same support but possibly different values of coefficients). Thus, when for every $i = 1, \dots, \ell$ we have $g(\mathbf{a}_i) \neq 0$, the method correctly computes all Gröbner bases.

Gröbner tracing fits naturally in Algorithm 6. Of course, we might not be able to ensure that for every $i = 1, \dots, \ell$ we have $g(\mathbf{a}_i) \neq 0$, but the probabilistic statement in Proposition 3.4 holds. In our implementation, $k$ is the field of integers modulo a prime, where the prime fits in 64 bits. We use the F4 algorithm [39] instead of the Buchberger's algorithm for more efficiency. Tracing in the F4 algorithm can be arranged similarly, with a difference that instead of S-polynomials we have rows of a Macaulay matrix. We use the package Groebner.jl [34], which implements the F4 algorithm and tracing as described above.

In Table 2, we illustrate the performance benefit of using Gröbner tracing to compute a single Gröbner basis in degrevlex over integers modulo a prime with our implementation. We report the running time of a classic F4 algorithm and of the F4 algorithm modified as described in the (Apply) bullet point (where we assume that the (Learn) part was already ran and we have access to its output). The running times are short so we compute them as averages of many runs. We observe a consistent speedup with the factor ranging from 1.8 to 6.1 across different problem sizes. This translates into a similar speedup for the whole Algorithm 6, whose running time is typically dominated by many runs of (Apply).

**Remark 5.3.** A different modification of the F4 algorithm that also exploits the low-rank structure of the matrices relies on probabilistic linear algebra over finite fields [97, Section 2]$^2$. However, unexpectedly, for many problems considered in this paper (in particular, from Table 2) our implementation of this method does not provide a significant speedup over the classic F4 algorithm.

**Remark 5.4.** A variant of the F4 algorithm that uses Gröbner tracing for computations over rational function fields has also been discussed in [72]. Their application is different, but we may compare the speedups obtained from tracing for a common example. For the Katsura-11 and Katsura-12 problems, over integers modulo a 32-bit prime, the quotients of the running times F4 algorithm F4 with tracing are 3.9 and 5.1, respectively, for the implementation in [72, Fig.3.]; for our implementation, they are 10.7 and 11.7, respectively. These appear to be of the similar order of magnitude; we also observe that tracing seems to be more beneficial for larger problems.

<footnote>$^2$We thank Roman Pearce for pointing out this observation to us.</footnote>
## 5.6 Comparison with other methods for computing Gröbner bases

While computing the full Gröbner basis with rational function coefficients is not the main focus of the present paper, Algorithm 6 allows doing this. Therefore, it may be interesting to compare this interpolation-based algorithm with other available methods.

<table><caption>Table 3: The running time of several different methods of computing Gröbner bases in $Q(x)[y]$ in degrevlex. Timeout means the computation was aborted after 24 hours. OOM means the computation was aborted after using more than 50 GB of RAM.</caption><thead><tr><th>Example</th><th>F4 (direct)</th><th>F4 (flat ring)</th><th>slimgb [19]</th><th>ffmodStd [17]</th><th>Algorithm 6</th></tr></thead><tbody><tr><td>Simson-3 [51]</td><td>1 s</td><td>1 s</td><td>0 s</td><td>0 s</td><td>0 s</td></tr><tr><td>Boku-I4 [16]</td><td>OOM</td><td>OOM</td><td>141 s</td><td>108 min</td><td>36 s</td></tr><tr><td>Param-1</td><td>OOM</td><td>3 s</td><td>3 s</td><td>20 min</td><td>OOM</td></tr><tr><td>Param-2</td><td>OOM</td><td>62 s</td><td>0 s</td><td>774 s</td><td>6 s</td></tr><tr><td>Moments-3 [59]</td><td>OOM</td><td>OOM</td><td>4 hours</td><td>Timeout</td><td>81 s</td></tr><tr><td>MAPK-5 [89]</td><td>14 s</td><td>OOM</td><td>3 s</td><td>Timeout</td><td>16 min</td></tr><tr><td>Pharm [32]</td><td>10 min</td><td>OOM</td><td>15 s</td><td>Timeout</td><td>100 s</td></tr><tr><td>Fujita [46]</td><td>0 s</td><td>0 s</td><td>0 s</td><td>52 s</td><td>0 s</td></tr><tr><td>Bilirubin [94]</td><td>1 s</td><td>5 s</td><td>0 s</td><td>113 s</td><td>0 s</td></tr><tr><td>EAIHRD [45]</td><td>OOM</td><td>OOM</td><td>Timeout</td><td>Timeout</td><td>Timeout</td></tr><tr><td>SLIQR [30]</td><td>0 s</td><td>OOM</td><td>0 s</td><td>127 s</td><td>0 s</td></tr><tr><td>LinComp2 [1]</td><td>OOM</td><td>OOM</td><td>OOM</td><td>Timeout</td><td>Timeout</td></tr></tbody></table>

We will be computing Gröbner bases of ideals in $Q(x)[y]$ in degrevlex ordering using the following methods:

1. The F4 algorithm applied directly in $Q(x)[y]$. We use the implementation in Groebner.jl [34], which is generic and has no special optimizations in this case.

2. The F4 algorithm applied in the flat ring $Q[x, y]$ with a suitable monomial order $<_{x,y}$. For $i = 1, 2$ for monomials $x^{\alpha_i}y^{\beta_i}$ with $\alpha_i \in \mathbb{N}_0^n$ and $\beta_i \in \mathbb{N}_0^t$, we have: $x^{\alpha_1}y^{\beta_1} <_{x,y} x^{\alpha_2}y^{\beta_2}$ if $y^{\beta_1} <_{d_{\text{rl}}} y^{\beta_2}$ or $y^{\beta_1} = y^{\beta_2}$ and $x^{\alpha_1} <_{d_{\text{rl}}} x^{\alpha_2}$, where $<_{d_{\text{rl}}}$ denotes the degrevlex monomial ordering. We use the implementation in Groebner.jl, which homogenizes the input generators in this case.

3. The slimgb algorithm [19] applied directly in $Q(x)[y]$. We use the implementation in Singular.

4. The ffmodStd algorithm [16, 17] implemented in Singular. This is an interpolation-based algorithm similar to Algorithm 6. ffmodStd reconstructs intermediate results to the rationals and performs interpolation over the rationals using the Cuyt-Lee's algorithm [29], while we carry out all computations (including interpolation) modulo a prime and reconstruct to the rationals only at the end. ffmodStd uses the built-in modStd command [66] for computing Gröbner bases in $Q[y]$ and does not use Gröbner tracing.

5. Algorithm 6, which we use to compute the full Gröbner basis.

In Table 3, we report the running time for each method for several examples. Each example is given by a list of generators of the ideal in $Q(x)[y]$. All implementations are run using only one CPU. Note that the methods 2., 4., and 5. are Monte-Carlo probabilistic. We made sure that the outputs of the implementations agree. Examples roughly fall into groups based on their provenance:

* Simson-3, distributed in the SymbolicData project [51]<fnref n="3" />. In the original slimgb publication, this system took the most time to compute among benchmark systems in $Q(x)[y]$ in degrevlex ordering. Now the computation is near instant with all methods.

* Boku-I4, from [16, Section 5.3.3]. This system appeared in the benchmarks of ffmodStd, where it took the longest time to run among benchmark systems in degrevlex. We modify the example by introducing a fictional variable in the ground ring (the system itself does not change) to disable the optimization that speeds up the computation over univariate rational function fields in ffmodStd<fnref n="4" />. In this example, the quotient of the running times $\frac{\text{ffmodStd}}{\text{Algorithm 6}}$ is around 180.

<footnote n="3">³We downloaded the systems' definitions from https://github.com/symbolicdata/data</footnote>

<footnote n="4">⁴With the optimization enabled, the running time of ffmodStd is 110 seconds</footnote>
* Param-1, Param-2. These were constructed from Katsura-5 by replacing some numerical coefficients with transcendental parameters. For Param-1, our implementation of Algorithm 6 has no chance to succeed, because the output contains a coefficient of total degree 46 in 3 variables, which is not possible to interpolate using our version of the Ben-Or and Tiwari algorithm modulo a 64-bit prime. For Param-2, the coefficients of the Gröbner basis contain large rationals (the numerators are 550 bits in size and the denominators are small), and our implementation uses 18 primes of size 64-bit to reconstruct the result.

* Moments-3, Example 7.11 from Section 7.3.

* Examples arising from the problem of structural identifiability (see Section 7.1). These have the structure of eOMS ideals (Notation 2.14). Tables 1 and 2 give an idea about the hardness of these examples, in particular we note that the Gröbner bases for these examples are typically almost linear and have very sparse coefficients.

Overall, computations tend to either finish quickly or not finish due to time and space constraints. As expected, straightforward direct methods (1. and 2.) require an unpredictable amount of memory and are thus difficult to use. The slimgb algorithm fares very well. Algorithm 6 is slower than slimgb for several examples, but is quicker for some others (e.g., Moments-3). Compared to ffmodStd, our implementation appears to be a hundred or more times more efficient; however, ffmodStd handles the Param-1 example, and our implementation does not. Finally, although none of the methods can compute the full Gröbner basis for the EAIHRD and LinComp2 examples, our Algorithm 6 readily tackles these examples using the approach from Section 5.4.

# 6 Experimental results

## 6.1 Assessing the quality of simplification

In this section, we use our implementation of the main algorithm (Algorithm 8) to simplify generators of subfields of rational function fields that arise in the problem of structural identifiability (see Section 7.1).

<table><caption>Table 4: Results of our simplification procedure: field generators before and after applying Algorithm 8</caption><thead><tr><th>Example</th><th>Original generators</th><th>Result of Algorithm 8</th></tr></thead><tbody><tr><td>Example 2.13</td><td>x<sub>1</sub><sup>2</sup> + x<sub>2</sub><sup>2</sup>, x<sub>1</sub><sup>3</sup> + x<sub>2</sub><sup>3</sup>, x<sub>1</sub><sup>4</sup> + x<sub>2</sub><sup>4</sup></td><td>x<sub>1</sub> + x<sub>2</sub>, x<sub>1</sub>x<sub>2</sub></td></tr><tr><td>genLV [110]</td><td>(β<sub>12</sub>r<sub>1</sub>r<sub>2</sub> - β<sub>22</sub>r<sub>1</sub><sup>2</sup>)/β<sub>12</sub>,<br/>(-β<sub>11</sub><sup>2</sup>β<sub>22</sub> + β<sub>11</sub>β<sub>12</sub>β<sub>21</sub>)/β<sub>12</sub>,<br/>(β<sub>12</sub>r<sub>2</sub> - 2β<sub>22</sub>r<sub>1</sub>)/β<sub>12</sub>, (β<sub>12</sub> + β<sub>22</sub>)/β<sub>12</sub>,<br/>(β<sub>11</sub>β<sub>12</sub>r<sub>2</sub> - 2β<sub>11</sub>β<sub>22</sub>r<sub>1</sub> + β<sub>12</sub>β<sub>21</sub>r<sub>1</sub>)/β<sub>12</sub>,<br/>(-β<sub>11</sub>β<sub>12</sub> + 2β<sub>11</sub>β<sub>22</sub> - β<sub>12</sub>β<sub>21</sub>)/β<sub>12</sub></td><td>r<sub>1</sub>, r<sub>2</sub>, β<sub>11</sub>, β<sub>21</sub>, β<sub>12</sub>/β<sub>22</sub></td></tr><tr><td>Pharm [32]</td><td>3533 functions, max. degree 31,<br/>2210 max. terms, size on disk 35 MB</td><td>n, k<sub>a</sub>, k<sub>c</sub>, b<sub>1</sub>, b<sub>2</sub>, a<sub>1</sub>, a<sub>2</sub></td></tr><tr><td>CGV1990 [24]</td><td>412 functions;<br/>indeterminates R and S appear<br/>in functions of degree at least 14</td><td>k<sub>3</sub>, k<sub>4</sub>, k<sub>6</sub>, k<sub>7</sub>, RV<sub>3</sub> + SV<sub>36</sub>,<br/>RV<sub>36</sub> + SV<sub>3</sub>/25,<br/>k<sub>5</sub>(RV<sub>36</sub> + SV<sub>36</sub>)/5,<br/>k<sub>5</sub>(V<sub>3</sub> + 5V<sub>36</sub>)/V<sub>3</sub></td></tr><tr><td>Fujita [46]</td><td>217 functions;<br/>indeterminates a<sub>1</sub>, a<sub>2</sub>, a<sub>3</sub> appear<br/>in functions of degrees 5, 3, 2, resp.</td><td>r<sub>22</sub>, r<sub>31</sub>, r<sub>41</sub>, r<sub>52</sub>, r<sub>61</sub>, r<sub>71</sub>,<br/>r<sub>81</sub>, r<sub>11</sub> - r<sub>12</sub> - r<sub>91</sub>, a<sub>1</sub>/r<sub>51</sub>,<br/>a<sub>3</sub>/r<sub>51</sub>, a<sub>2</sub>/r<sub>51</sub>, r<sub>21</sub>/r<sub>51</sub></td></tr><tr><td>JAK-STAT [109]</td><td>1300 functions;<br/>indeterminates t<sub>15</sub> and t<sub>21</sub> appear<br/>in functions of degree at least 4</td><td>t<sub>1</sub>, t<sub>2</sub>, t<sub>3</sub>, t<sub>4</sub>, t<sub>5</sub>, t<sub>6</sub>, t<sub>7</sub>, t<sub>8</sub>, t<sub>9</sub>,<br/>t<sub>10</sub>, t<sub>12</sub>, t<sub>13</sub>, t<sub>14</sub>, t<sub>16</sub>, t<sub>18</sub>, t<sub>19</sub>,<br/>t<sub>20</sub>, t<sub>17</sub>t<sub>22</sub>, t<sub>15</sub>t<sub>21</sub>, t<sub>11</sub>t<sub>21</sub></td></tr></tbody></table>

Our benchmark suite consists of 53 examples, which correspond to problems found in literature. Many of these examples have been collected by Barreiro and Villaverde in [8, 111] and in the works referenced therein. For each example, we construct a list of rational functions, the generating set to be simplified, and pass it to our main algorithm. We run the algorithm with the same parameters for all examples: Gröbner bases are computed in degrevlex ordering, Algorithm 7 uses the degree bound $\delta := 3$, and (Step 6) orders the generators using the strategy from Remark 3.2 and does not enforce minimality.
In Table 4, we list results for a selection of examples (the rest are given in Appendix A) and also include Example 2.13 to clarify the format. In this section, for brevity, the degree of a rational function means the sum of the degrees of its numerator and denominator. Since we did not state the simplification problem in a fully formal way, we have no single qualitative measure to assess the quality of simplification. Therefore, we use some natural related measures (such as the number of generators, their degrees, and further criteria from Remark 3.2) in combination with manual inspection. Analyzing the overall results, we make the following observations, which we will illustrate using the examples from the table.

* Algorithm 8 drastically reduces the size of the generating sets. Note that this does not amount to simply removing excessive generators: we observe that even the *minimal* degrees of generators for some of the variables may be reduced significantly.

* In particular, we see instances of nontrivial rewriting. In the genLV example, the original generating set has a common denominator $\beta_{12}$, and the variables $r_1, r_2$ always appear together and in fractions with the numerator of degree at least two. In the simplified generating set there is only one fraction, and $r_1$ and $r_2$ are separated. In the Pharm example, every variable appears in the original generating set in an element of degree at least two, but the result of simplification is of degree one.

* Simplification reveals that some subfields can be generated by polynomials. In the JAK-STAT example, the original generating set contains rational functions with denominators of degree up to six, but the subfield can be generated by simple polynomials only. The sorting criteria from Remark 3.2 facilitate such simplifications by favoring elements with small denominators. Overall, in 26 out of 53 examples, the original generating set contained rational functions with non-constant denominators, but the simplification revealed that the subfield can be generated by polynomials only.

* In every but one example (Bilirubin), the resulting generating set is algebraically independent over $\mathbb{C}$. Independence implies that the generating set contains the fewest elements possible [84, Ch. VIII, Theorem 1.1.].

Each example was run with a 1-hour timeout and a 50 GB RAM limit. The median running time was 20 seconds. The longest running example was Ovarian_follicle with the running time of 23 minutes. One example, Covid-2, did not finish getting stuck in (Step 2) and (Step 3) of Algorithm 8. This example stands apart because the coefficients of the OMS ideal are considerably larger than the input, making interpolation infeasible. Even in the case of succesfull computation of these coefficients, it is not clear whether they would be useful for simplification. More precisely, the largest degree in the input is 29, while the coefficients of the Gröbner basis of the OMS ideal that generate the subfield are of degree 229. This is not typical for examples coming from this domain, and we do not what makes this example special.

From this experiment, we conclude that our method effectively deals with most of the real-world examples in the benchmark suite. In Section 7, we will discuss the usefulness of such simplifications in applications, including but not limited to structural identifiability for ODEs. Additionally, given that we handle most of the benchmark problems, a natural direction would be to apply our algorithm to harder examples from the same area. One current limitation is that producing generating sets for such examples using the approach from [37] (see Section 7.1) is itself computationally demanding.

## 6.2 Comparison with other methods for simplification

We will compare our method with the simplification algorithm from [105] implemented in Maple<fnref n="5" />. This method also uses the coefficients of OMS ideals, but computes them directly instead of using interpolation. To pick the simplest elements from a generating set, it orders the functions by the length of their string representation. It does not implement Algorithm 7. This method has been developed with the applications in mind that are similar to ours.

<footnote n="5">⁵We used the implementation from https://github.com/pogudingleb/AllIdentifiableFunctions</footnote>
<table><thead><tr><th>Example</th><th>Result of [105]</th><th>Result of Algorithm 8</th></tr></thead><tbody><tr><td>Biohydrogenation [96]</td><td>k<sub>5</sub>, k<sub>6</sub>, k<sub>7</sub>, k<sub>0</sub><sup>2</sup>, k<sub>10</sub>/k<sub>9</sub>,<br>(k<sub>10</sub>k<sub>9</sub> + 2k<sub>8</sub>k<sub>9</sub>)/k<sub>10</sub></td><td>k<sub>5</sub>, k<sub>6</sub>, k<sub>7</sub>, k<sub>0</sub><sup>2</sup>, k<sub>9</sub>k<sub>10</sub>,<br>2k<sub>8</sub> + k<sub>10</sub></td></tr><tr><td>SEIR1 [131]</td><td>γ, β/ψ, γψ - ψ - v,<br>β<sup>2</sup>(γ - v)/(γψ - ψ)</td><td>γ, β/ψ, ψv - ψ - v,<br>γψ - ψv</td></tr><tr><td>SLIQR [30]</td><td>σ, β, N, α + γ, (ησ + α - σ)/αη,<br>η<sup>2</sup>γσ<sup>2</sup> + α<sup>2</sup>ησ + 3αηγσ - αησ<sup>2</sup> - 2ηγσ<sup>2</sup> + α<sup>2</sup>γ - 2αγσ + γσ<sup>2</sup>/ησ + α - σ</td><td>σ, β, N, α + γ, αηγ,<br>αγ + ηγσ - γσ</td></tr></tbody></table>

<table><caption>Table 5: The differences of results obtained by the method from [105] and our implementation of Algorithm 8 on some cherry-picked examples</caption><tbody><tr><td rowspan="2">Bilirubin [94]</td><td>k<sub>01</sub>, k<sub>12</sub> + k<sub>13</sub> + k<sub>14</sub>, k<sub>21</sub> + k<sub>31</sub> + k<sub>41</sub>,<br>k<sub>12</sub><sup>2</sup>k<sub>13</sub> - k<sub>12</sub><sup>2</sup>k<sub>14</sub> - k<sub>12</sub>k<sub>13</sub><sup>2</sup> + k<sub>12</sub>k<sub>14</sub><sup>2</sup> + k<sub>13</sub><sup>2</sup>k<sub>14</sub> - k<sub>13</sub>k<sub>14</sub><sup>2</sup>,<br>k<sub>12</sub>k<sub>31</sub> - k<sub>12</sub>k<sub>41</sub> - k<sub>13</sub>k<sub>21</sub> + k<sub>13</sub>k<sub>41</sub> + k<sub>14</sub>k<sub>21</sub> - k<sub>14</sub>k<sub>31</sub>,</td><td>k<sub>01</sub>, k<sub>12</sub>k<sub>13</sub>k<sub>14</sub>, k<sub>21</sub>k<sub>31</sub>k<sub>41</sub>,<br>k<sub>12</sub> + k<sub>13</sub> + k<sub>14</sub>, k<sub>21</sub> + k<sub>31</sub> + k<sub>41</sub>,<br>k<sub>21</sub>k<sub>31</sub> + k<sub>21</sub>k<sub>41</sub> + k<sub>31</sub>k<sub>41</sub>,</td></tr><tr><td>k<sub>21</sub><sup>2</sup>k<sub>31</sub> - k<sub>21</sub><sup>2</sup>k<sub>41</sub> - k<sub>21</sub>k<sub>31</sub><sup>2</sup> + k<sub>21</sub>k<sub>41</sub><sup>2</sup> + k<sub>31</sub><sup>2</sup>k<sub>41</sub> - k<sub>31</sub>k<sub>41</sub><sup>2</sup>,<br>k<sub>12</sub>k<sub>31</sub> - k<sub>12</sub>k<sub>41</sub> - k<sub>13</sub>k<sub>21</sub> + k<sub>13</sub>k<sub>41</sub> + k<sub>14</sub>k<sub>21</sub> - k<sub>14</sub>k<sub>31</sub>,</td><td>k<sub>12</sub>k<sub>13</sub> + k<sub>12</sub>k<sub>14</sub> + k<sub>13</sub>k<sub>14</sub>,<br>k<sub>12</sub>k<sub>31</sub> + k<sub>12</sub>k<sub>41</sub> + k<sub>13</sub>k<sub>21</sub> + k<sub>13</sub>k<sub>41</sub>,<br>+k<sub>14</sub>k<sub>21</sub> + k<sub>14</sub>k<sub>31</sub>, k<sub>21</sub> + k<sub>31</sub> + k<sub>41</sub></td></tr><tr><td>-k<sub>12</sub>k<sub>31</sub><sup>2</sup> + k<sub>12</sub>k<sub>41</sub><sup>2</sup> + k<sub>13</sub>k<sub>21</sub><sup>2</sup> - k<sub>13</sub>k<sub>41</sub><sup>2</sup> - k<sub>14</sub>k<sub>21</sub><sup>2</sup> + k<sub>14</sub>k<sub>31</sub><sup>2</sup>,<br>k<sub>12</sub>k<sub>31</sub> - k<sub>12</sub>k<sub>41</sub> - k<sub>13</sub>k<sub>21</sub> + k<sub>13</sub>k<sub>41</sub> + k<sub>14</sub>k<sub>21</sub> - k<sub>14</sub>k<sub>31</sub><sup>2</sup>,</td><td></td></tr><tr><td>-k<sub>2</sub><sup>2</sup>k<sub>31</sub> + k<sub>2</sub><sup>2</sup>k<sub>41</sub> + k<sub>2</sub><sup>2</sup>k<sub>31</sub><sup>2</sup> - k<sub>2</sub><sup>2</sup>k<sub>41</sub> - k<sub>2</sub><sup>2</sup>k<sub>41</sub><sup>2</sup> + k<sub>2</sub><sup>2</sup>k<sub>41</sub><sup>2</sup> + k<sub>2</sub><sup>2</sup>k<sub>31</sub><sup>2</sup></td><td></td></tr><tr><td>k<sub>12</sub>k<sub>31</sub> - k<sub>12</sub>k<sub>41</sub> - k<sub>13</sub>k<sub>21</sub> + k<sub>13</sub>k<sub>41</sub> + k<sub>14</sub>k<sub>21</sub> - k<sub>14</sub>k<sub>31</sub><sup>2</sup></td><td></td></tr></tbody></table>

We run the method from [105] on our benchmark suite, with the same time and memory constraints as in Section 6.1. Overall, this method does not finish for 8 examples$^6$, obtains a result similar to ours for 36 examples, and in the remaining 9 examples our method’s results seem better. To corroborate the latter claim, we present the differences for some of the results in Table 5. For Biohydrogenation, SEIR1, and SLIQR, the methods produce generating sets with the same number of elements; however, our method produces elements of smaller degrees. In the Bilirubin example, our output is visually simpler but contains more elements (8 vs. 7); this also means that the result of [105] is algebraically independent, while ours is not. Of course, in some use cases, the result of [105] may be preferable.

We conclude that our method appears to be an improvement over [105] both in terms of better performance and for producing simpler generating sets.

**Remark 6.1 (Computing the algebraic closure).** In some applications, such as structural parameter identifiability, it may be of interest to consider not only a subfield $\mathcal{E}$ of $\mathbb{Q}(x)$ but also its relative algebraic closure $\mathcal{E}^{\text{rel}}$ in $\mathbb{Q}(x)$. An approach for finding a simple transcendence basis for $\mathcal{E}^{\text{rel}}$ based on a construction reminiscent to OMS ideals was proposed in a series of papers [92, 93, 94] and implemented in the COMBOS web-application [94]. A completely different approach for finding simple elements in $\mathcal{E}^{\text{rel}}$, based on linear algebra and avoiding Gröbner basis computations, was proposed in [91, Section 5]. A direct comparison of these methods to ours is complicated: for example, all the fields in Table 5 have the full $\mathbb{Q}(x)$ as the algebraic closure, so $x$ itself could be taken as a simple transcendence basis of $\mathcal{E}^{\text{rel}}$.

**Example 6.2.** Some simplified generators appearing in the literature can be sometimes simplified even further using our algorithm. Consider, for example, the following field generators produced by COMBOS in [94, Test-Model 2, $n = 3$]:

$$ \begin{aligned} & p_{11}, p_{21}, p_{22} + p_{33}, p_{23}/p_{13}, p_{22} - p_{12}p_{23}/p_{13}, \\ & p_{12}p_{23}p_{33}/p_{13} - p_{23}p_{32} + p_{12}p_{22}p_{23}/p_{13} - p_{22}^{2}, \\ & p_{11}p_{23}p_{33}/p_{13} - p_{23}p_{31} + p_{11}p_{22}p_{23}/p_{13} - p_{21}p_{22}. \end{aligned} $$

We can either feed these generators to our simplification algorithm, or use the approach from Section 7.1. Both produce the same result, which is visually simpler and contains elements of lower degrees:

$$ p_{11}, p_{21}, p_{22} + p_{33}, p_{13}/p_{23}, p_{22}p_{33} - p_{23}p_{32}, p_{21}p_{33} - p_{23}p_{31}, p_{12}p_{21} + p_{13}p_{31}. $$

# 7 Applications

In this section, we present three sets of examples that showcase the simplifications achieved by our algorithm. We start with examples arising in the study of structural identifiability of ODE models in Section 7.1. This application domain has been our original motivation for this work. Then we show that our algorithm

<footnote>$^6$The examples where the method does not finish also do not finish with the time limit increased ten-fold</footnote>
produces interesting results applied in other contexts such as identifiability and observability of discrete-time models (Section 7.2) and rational invariant theory (Section 7.3).

In addition to the situations where one is interested in simpler expressions for generators (such as the case studies below), we expect that having more concise generators can improve the performance of other algorithms on rational function fields (such as membership testing, see also Remark 4.6).

## 7.1 Structural identifiability of ODE models

The main motivation for the algorithmic developments presented in this paper comes from the problem of computing identifiable functions of dynamical models. We start the section with a high-level explanation of the way generator simplification question appears in modeling and its importance. For further technical details, we refer to [92, 105, 106]. Then we will showcase simplifications achieved by our algorithm using several systems from modeling literature and demonstrate that the simplified generating set is better suited to derive insights about the structure of the model. However, since mathematical modeling is not a central topic in this paper, we will not get deep into the domain-specific details.

We consider a parametric system of differential equations

$$x' = g(\mu, x), \qquad (8)$$

where $x$ is a vector of unknown functions of time, $\mu$ is a vector of scalar parameters, and $g$ is a vector of rational functions from $\mathbb{C}(x, \mu)$ defining the dynamics. The system may also involve external inputs, we will omit them to simplify the presentation. In a typical experimental setup, it is not feasible to measure the values of all the coordinates of $x$. Thus, data is available only for some functions of $x$ and $\mu$ typically referred to as outputs of the model. For the sake of simplicity, we will consider the case of a single output $y = f(\mu, x)$.

A parameter $\mu_i$ or, more generally, a function $h(\mu)$ of parameters is (structurally) identifiable if its value can be uniquely reconstructed from the output data assuming continuous and noise-free measurements (see [57, Section 2] for more details). Since identifiability is a prerequisite for meaningful parameter estimation, it is an important and desirable property. The set of identifiable functions can be shown to be closed under addition, multiplication, and division, so they form the field of identifiable functions. The elements of this field (also called identifiable combinations) are exactly the quantities whose values can be in principle inferred from data, so computing a convenient representation of this field is a natural question.

**Example 7.1.** Consider a model

$$x_1' = \mu_1 x_2, \quad x_2' = -\mu_2 x_1,$$

and assume that the observation is $y = x_1$, the first coordinate. First, we observe that the scaling $x_2 \to \lambda x_2, \mu_2 \to \lambda \mu_2, \mu_1 \to \frac{1}{\lambda} \mu_1$ preserves the output function for any nonzero constant $\lambda$. Thus, neither of $\mu_1$ and $\mu_2$ can be identified even from perfect knowledge of $y$. On the other hand, $y'' = -\mu_1 \mu_2 y$, so the product $\mu_1 \mu_2$ can be found as long as $y$ is not identically zero.

The classical approach for finding a generating set of the field of identifiable functions going back to [104] is the following. First, we compute the minimal polynomial differential equation (first with respect to the order, then degree) with coefficients in $\mathbb{C}(\mu)$ satisfied by $y$ (called the input-output equation). By normalizing one of the coefficients to one, we can write it as

$$m_0(y, \dots, y^{(h)}) + \sum_{i=1}^{M} c_i(\mu) m_i(y, \dots, y^{(h)}) = 0, \qquad (9)$$

where $m_0, \dots, m_M$ are distinct monomials in $y, y', \dots, y^{(h)}$ and $c_1, \dots, c_M \in \mathbb{C}(\mu)$. Then, under certain mild assumptions, the functions $c_1, \dots, c_M$ generate the whole field of identifiable functions [106, Corollary 1].

**Example 7.2 (Example 7.1 continued).** The minimal differential equation satisfied by $y$ in this example is $y'' + \mu_1 \mu_2 y = 0$. Therefore, [106, Corollary 1] implies that identifiable functions for this model are precisely rational functions in $\mu_1 \mu_2$.

Unlike the toy example above, in more realistic models, the functions $c_1, \dots, c_M$ are often large and complicated (see Appendix A); this makes such generating set practically useless. Thus, the challenge is, given a generating set $c_1, \dots, c_M$, to produce a simpler one which can later be understood and interpreted by a modeler. Below we will demonstrate our simplification algorithm on several fields coming from the structural identifiability problem.
**Example 7.3 (Compartmental model).** Consider the following linear model (which was employed, for example, for modeling bilirubin kinetics, see [11, Figure 1] and [94, Example 4]):

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

In this model, we have

$$x = [x_1, x_2, x_3] \quad \text{and} \quad \mu = [\mu_{01}, \mu_{12}, \mu_{13}, \mu_{21}, \mu_{31}],$$

observations are taken to be $y(t) = x_1(t)$, and $u(t)$ is an external input. Due to the linearity of the model, there will be not too many coefficients in the input-output equation (9), and taking them directly would yield the following generators of the field of identifiable functions:

$$\mu_{12}+\mu_{13}, \mu_{12}\mu_{13}, \mu_{01}\mu_{12}\mu_{13}, \mu_{01}+\mu_{12}+\mu_{13}+\mu_{21}+\mu_{31}, \mu_{01}\mu_{12}+\mu_{01}\mu_{13}+\mu_{12}\mu_{13}+\mu_{12}\mu_{31}+\mu_{13}\mu_{21}.$$

The simplified generating set produced by our Algorithm 8 is

$$\mu_{01}, \mu_{21} + \mu_{31}, \mu_{21}\mu_{31}, \mu_{12} + \mu_{13}, \mu_{12}\mu_{13}, \mu_{12}\mu_{31} + \mu_{13}\mu_{21}. \quad (11)$$

This generating set reveals several important properties of the model:

* identifiability of $\mu_{01}$ is made explicit (in the original generating set, $\mu_{01}$ could be obtained as the quotient of the second and third generators);

* the symmetry with respect to the exchange of the second and third compartments ($x_2$ and $x_3$) also becomes explicit; it is, furthermore, clear that the field is precisely the field of invariants with respect to the simultaneous exchange ($\mu_{12}, \mu_{21} \leftrightarrow (\mu_{13}, \mu_{31})$).

We note, however, that the simplified generating set contains one more element than the original one (6 vs. 5), so the generators are not independent. The minimality and independence can be restored by applying Algorithm 5, the resulting set will be (11) without $\mu_{12}\mu_{13}$. While this set is smaller, it does not reveal the symmetry $x_2 \leftrightarrow x_3$ as explicitly, so, depending on the goals, one of these two outputs may be preferable.

**Example 7.4 (SIRT epidemiological model).** We consider an epidemiological model describing a spread of infection which takes into account that the individuals receiving treatment are less infectious [123, Eq. (2.3)]:

$$\left\{ \begin{array}{l} S'(t) = -\frac{\beta S(t) I(t)}{N} - \frac{\delta \beta S(t) T(t)}{N}, \\ I'(t) = \frac{\beta S(t) I(t)}{N} + \frac{\delta \beta S(t) T(t)}{N} - (\alpha + \gamma) I(t), \\ T'(t) = \gamma I(t) - \nu T(t). \end{array} \right. \quad \text{where} \quad \begin{array}{l} \mathbf{x} = [S, I, T] \\ \mu = [\alpha, \beta, \gamma, \delta, \nu, N] \end{array} \quad (12)$$

In this model, $S$ and $I$ denote the number of susceptible and infectious individuals, respectively, while $T$ stands for the number of individuals receiving treatment. The observation is taken to be $y(t) = T(t)$.

The generating set of the field of identifiable functions extracted directly from the input-output equation (9) consists of twelve rational elements: polynomials up to degree six and rational functions with linear denominator and numerator up to degree five. Algorithm 8 produces the following generating set for the same field:

$$\alpha + \gamma + \nu, \quad \nu(\alpha + \gamma), \quad \nu + \delta\gamma, \quad \frac{N\gamma}{\beta}. \quad (13)$$

The first two generators are the elementary symmetric functions in $\alpha + \gamma$ and $\nu$ hinting at a symmetry between the outflow rates of the $I$ and $T$ compartments of the model. Indeed, such a symmetry is a recurring motif (see, e.g., [56, A.3, Eq. (9)] and [107, Example 4.1]).

To facilitate the interpretation of the third generator in (13), we observe that the nonlinear terms in the first two equations (12) can be factored as $\frac{\beta S(t)}{N}(I(t) + \delta T(t))$ suggesting $\tilde{I}(t) := I(t) + \delta T(t)$ ("total
infectiousness”) as a natural quantity to consider. Rewriting the equation for $T$ in terms of $T$ and $\tilde{I}(t)$, we get

$$T'(t) = \gamma \tilde{I}(t) - (\nu + \gamma \delta)T(t),$$

where the third generator from (13) appears as an outflow rate.

**Example 7.5 (SLIQR epidemiological model).** Consider an epidemiological model describing a spread of infection which takes into account the latent period and possible remissions ([31, Equation (4)]):

$$ \left\{ \begin{aligned} S'(t) &= -\frac{\beta S(t) I(t)}{N} - v(t) \frac{S(t)}{N}, \\ L'(t) &= \frac{\beta S(t) \tilde{I}(t)}{N} - \alpha L(t), \\ I'(t) &= \alpha L(t) + \sigma Q(t) - \gamma I(t), \\ Q'(t) &= \gamma(1 - \eta)I(t) - \sigma Q(t) \end{aligned} \right. \quad \text{where} \quad \begin{array}{l} \mathbf{x} = [S, L, I, Q] \\ \mathbf{\mu} = [\alpha, \beta, \gamma, \eta, \sigma, N] \end{array} \qquad (14) $$

In this model, $S$ and $I$ denote the number of susceptible and infectious individuals, respectively. $L$ and $Q$ stand for the number of individuals in the latent and remission stages, respectively. Furthermore, $v(t)$ is an external input (the proportion of vaccinated individuals) and the observation is taken to be the prevalence (the proportion of infectious individuals), that is, $y(t) = \frac{I(t)}{N}$.

The generating set of the field of identifiable functions extracted directly from the input-output equation (9) consists of 40 polynomials with the degree reaching eight. Algorithm 8 produces the following set of generators for the same field:

$$\beta, \sigma, N, \alpha + \gamma, \alpha\eta\gamma, \alpha\gamma + (\eta - 1)\gamma\sigma.$$

First, this generating set explicitly features identifiable parameters $\beta, \sigma$, and $N$. Two of the remaining generators appear naturally if we introduce a new quantity, $C(t) := I'(t)$, the number of new infected cases, a quantity frequently used in epidemiology (also called incidence). A direct computation shows that the differential equation for $C(t)$ will be

$$C'(t) = \frac{\alpha\beta S(t)I(t)}{N} - (\alpha + \gamma)C(t) - (\alpha\gamma + (\eta - 1)\gamma\sigma)I(t) + \sigma(\alpha - \sigma)Q(t).$$

This equation explicitly involves the forth and sixth generators as the flow rates.

## 7.2 Identifiability and observability for discrete-time systems

In this section we will apply our algorithm for the problems of similar flavor as in the previous subsection but for models in discrete time. The general setup is very similar although the computation will be quite different. Specifically, we consider a parametric system of difference equations (cf. (8)):

$$\mathbf{x}(t + 1) = \mathbf{g}(\mathbf{\mu}, \mathbf{x}(t)), \qquad (15)$$

where $\mathbf{x}(t)$ is a vector of time-series (for $t = 0, 1, 2, \dots$), $\mathbf{\mu}$ is a vector of unknown parameters, and $\mathbf{g}$ is a vector of rational functions from $\mathbb{C}(\mathbf{x}, \mathbf{\mu})$. We assume that we observe an output time series $y(t) = f(\mathbf{\mu}, \mathbf{x}(t))$, and are interested in identifiable functions in this context. In other words, we are given a sequence $y(0), y(1), y(2), \dots$ or a finite part of it, and we want to find out which functions of parameters $\mathbf{\mu}$ and initial conditions $\mathbf{x}(0)$ we can reconstruct. For generic values of $\mathbf{\mu}$ and $\mathbf{x}(0)$ (so that none of the relevant denominators vanish), these functions form a field which we will call the field of identifiable functions (in the literature, for initial conditions $\mathbf{x}(0)$, term *observable functions* is also often employed).

Despite of the similarity of the general question to the one in Section 7.1, the approach via input-output equations does not work in the difference case as well as it does in the differential (see [88] for a discussion). Therefore, we will use the following more direct method (cf. [76]). We define rational functions $\mathbf{g}^{\circ j}(\mathbf{\mu}, \mathbf{x}(0))$ by

$$\mathbf{g}^{\circ 0}(\mathbf{\mu}, \mathbf{x}(0)) := \mathbf{x}(0) \quad \text{and} \quad \mathbf{g}^{\circ (j+1)}(\mathbf{\mu}, \mathbf{x}(0)) := \mathbf{g}(\mathbf{\mu}, \mathbf{g}^{\circ j}(\mathbf{\mu}, \mathbf{x}(0))) \text{ for } j \ge 0.$$

Comparing with (15), we observe that $\mathbf{x}(j) = \mathbf{g}^{(j)}(\mathbf{\mu}, \mathbf{x}(0))$. Thus, $y(j)$ can be written as a rational function $f(\mathbf{\mu}, \mathbf{g}^{(j)}(\mathbf{\mu}, \mathbf{x}(0)))$, and the field of functions identifiable from $y(0), \dots, y(\ell)$ is generated by these rational functions for $j = 0, \dots, \ell$. However, as we will see below, these rational functions become complicated very quickly, so simplification is a key step in finding an interpretable set of generators.
**Example 7.6 (General fractional-linear model).** We will consider a one-dimensional model with dynamics defined by a general fractional-linear function:

$$x(t + 1) = \frac{\mu_1 x(t) + \mu_2}{\mu_3 x(t) + \mu_4}, \quad y(t) = x(t).$$

We start with computing simple generators for the field generated by the values of the output at $t = 0, 1, 2, 3$. These values can be written as rational functions in $\mu$, $x(0)$ as follows:

$$y(0) = x(0), \quad y(1) = \frac{\mu_1 x(0) + \mu_2}{\mu_3 x(0) + \mu_4}, \quad y(2) = \frac{(\mu_1^2 + \mu_2\mu_3)x(0) + \mu_1\mu_2 + \mu_2\mu_4}{(\mu_1\mu_3 + \mu_3\mu_4)x(0) + \mu_2\mu_3 + \mu_4^2}$$

$$y(3) = \frac{(\mu_1^3 + 2\mu_1\mu_2\mu_3 + \mu_2\mu_3\mu_4)x(0) + \mu_1^2\mu_2 + \mu_1\mu_2\mu_4 + \mu_2^2\mu_3 + \mu_2\mu_4^2}{(\mu_1^2\mu_3 + \mu_1\mu_3\mu_4 + \mu_2\mu_3^2 + \mu_3\mu_4^2)x(0) + \mu_1\mu_2\mu_3 + 2\mu_2\mu_3\mu_4 + \mu_4^3}$$

Our simplification algorithm finds the following much simpler generators for the same field:

$$x(0), \frac{\mu_1}{\mu_4}, \frac{\mu_2}{\mu_4}, \frac{\mu_3}{\mu_4}.$$

One can show that taking further time points will not add any new information since the fraction in the right-hand side is invariant under multiplying all the $\mu$'s by the same number. On the other hand, the output of our algorithm suggests that the fields obtained by taking fewer time points are strictly smaller. Therefore, the field of identifiable functions is generated by $x(0), \frac{\mu_1}{\mu_4}, \frac{\mu_2}{\mu_4}, \frac{\mu_3}{\mu_4}$ and four time points are enough to find these values.

**Example 7.7 (Logistic growth).** The following logistic growth model is used, for example, in population modeling [3, p. 12]:

$$x(t + 1) = x(t) \left( 1 + \mu_1 \left( 1 - \frac{x(t)}{\mu_2} \right) \right).$$

We will further assume that what is observed is an unknown but fixed fraction of the population, that is, $y(t) = \mu_3 x(t)$. In order to find identifiable functions, we will consider the field generated by $y(0), y(1)$, and $y(2)$, where the latter is of the form $\frac{P(x(0), \mu)}{\mu_2^3}$ with $\deg P = 8$. In contrast to this, our algorithm finds the following generators for the same field

$$\mu_3 x(0), \mu_1, \mu_2 \mu_3.$$

By repeating this computation with $y(4)$, one can show that adding new time points does not change the generators. Since the field is smaller if only one or two time points are taken, we conclude that the three functions above generate the field of identifiable functions and knowledge of only three data points is generically enough to reconstruct them.

**Example 7.8 (SIS/SIR epidemiological model).** The following model combines the dynamics of discrete-time SIS and SIR models [2, Sections 3 and 4]:

$$\left\{
\begin{array}{l}
S(t + 1) = S(t)(1 - \beta I(t)) + \alpha I(t), \\
I(t + 1) = I(t)(1 - \alpha - \gamma + \beta S(t)), \\
R(t + 1) = R(t) + \gamma I(t), \\
y(t) = \eta I(t),
\end{array}
\right.
\quad \text{where} \quad
\begin{array}{l}
\mathbf{x} = [S, I, R] \\
\mathbf{\mu} = [\alpha, \beta, \gamma, \eta]
\end{array}
\qquad (16)$$

As in the previous examples, we compute the field generated by the values $y(0), y(1), \dots$ and observe that after $y(0), \dots, y(3)$ the field stabilizes. Simplifying $y(0), \dots, y(3)$, where the latter is a polynomial of degree 16 with 102 terms, we obtain a much simpler generating set of identifiable functions:

$$\gamma, \eta I(0), \beta I(0), \beta S(0) - \alpha.$$

These simpler expressions also admit biological interpretations: $\eta I(0)$ is the observed number of infections at $t = 0$, $\beta I(0)$ is essentially the total infection rate at $t = 0$, and $\alpha - \beta S(0)$ is the total rate of the flow from the I-compartment to the S-compartment.

## 7.3 Fields of invariants

Let $G$ be a group and consider its action on an $n$-dimensional $k$-vector space $V$ (in other words, a homomorphism $G \to GL(V)$). This action naturally defines the action of $G$ on the ring of polynomials $k[V]$
and the field of rational fractions $k(V)$. Invariant theory studies the invariants of this action: the ring of polynomial invariants $k[V]^G$ and the field of rational invariants $k(V)^G$. Historically, the invariant theory was focused on the polynomial invariants (see, e.g. [128]), but in the recent years rational invariants have attracted attention (see, e.g. [7, 15, 18, 61, 62, 63, 81]). From the applications point of view, invariants are useful to characterize data (signals or images) which is defined up to a group action [7]. The degrees of the generators are related to the sampling complexity of the reconstruction problem for such data [7], and the generators themselves can be used to distinguish different orbits with respect to the action [42].

The first algorithm for computing a generating set of $k(V)^G$ was proposed by Müller-Quade and Beth in [99]. It was further extended by Hubert and Kogan in [64], in particular, by allowing to provide a cross-section [64, Section 3.1] of the group action as an extra input (for some other generalizations, see [78]). Providing a cross-section typically reduces the computation time and results in a simpler generating set. In the following examples we show that the methods and software described in this paper may simplify further the generating sets produced by the state of the art algorithms.

For the computations below, we used an implementation of the Hubert-Kogan algorithm from the AIDA [60] Maple package.

**Example 7.9 (S<sub>3</sub>-action).** Consider the action of the symmetric group $S_3$ on 6-dimensional space (inspired by [50, Section 4.1]) such that a permutation $\sigma \in S_3$ acts on a vector $(a_1, a_2, a_3, b_1, b_2, b_3) \in \mathbb{R}^6$ as follows:

$$ \sigma(a_1, a_2, a_3, b_1, b_2, b_3) = (a_{\sigma(1)}, a_{\sigma(2)}, a_{\sigma(3)}, (-1)^{\sigma}b_{\sigma(1)}, (-1)^{\sigma}b_{\sigma(2)}, (-1)^{\sigma}b_{\sigma(3)}), $$

where $(-1)^\sigma$ denotes the sign of the permutation. Finite groups do not admit cross-sections, so we run the Hubert-Kogan algorithm without providing any additional information. The resulting generating set contains a natural symmetric polynomial $a_1 + a_2 + a_3$ and 33 rational functions with total degrees ranging from five to seven. We apply our simplification algorithm to this set and obtain a simpler generating set consisting of seven polynomials only:

$$ \begin{aligned} & a_1 + a_2 + a_3, \quad a_1 a_2 + a_1 a_3 + a_2 a_3, \quad a_1 a_2 a_3, \quad b_1 b_2 + b_1 b_3 + b_2 b_3, \\ & b_1^2 + b_2^2 + b_3^2, \quad a_1 b_2 - a_1 b_3 - a_2 b_1 + a_2 b_3 + a_3 b_1 - a_3 b_2, \quad a_1 b_2 b_3 + a_2 b_1 b_3 + a_3 b_1 b_2. \end{aligned} $$

**Example 7.10 (Triangular action with swap).** Consider a space $V$ of at most quadratic bivariate polynomials without constant term $a_{01}x + a_{10}y + a_{02}x^2 + 2a_{11}xy + a_{20}y^2$. Group $T_2$ of nonsingular upper-triangular $2 \times 2$ matrices acts on $V$ via its action on the $(x, y)$-space. We define the action of $T_2 \times (\mathbb{Z}/2\mathbb{Z})$ on $V \oplus V \cong \mathbb{R}^{10}$ by $T_2$ acting on the summands and $(\mathbb{Z}/2\mathbb{Z})$ swapping the two copies of $V$. The coordinates on the second copy of $V$ will be denoted by $b_{01}, b_{10}, b_{02}, b_{11}, b_{20}$.

We apply the Hubert-Kogan algorithm with the cross-section $a_{11} = 0, a_{02} = 1$. This cross-section is chosen as the one with the simplest output among all cross-sections consisting of the equations of the form $x = 0$ or $x = 1$, where $x$ is a variable$^7$, for which the computation finishes within 10 minutes. The resulting generating set contains $\frac{a_{11}b_{02}-a_{02}b_{11}}{a_{02}-b_{02}}$ and $\frac{a_{02}^2+b_{02}^2}{a_{02}b_{02}}$ as well as 18 rational function of total degrees ranging from 5 to 11. Our simplification algorithm produces the following generating set consisting of five polynomials and four rational functions (including the two simple elements of the original generating set) with degrees ranging from two to four:

$$ \begin{aligned} & a_{20}b_{02} - 2a_{11}b_{11} + a_{02}b_{20}, \quad a_{20}a_{02} - a_{11}^2 + b_{20}b_{02} - b_{11}^2, \\ & a_{10}^2a_{02} - 2a_{10}a_{01}a_{11} + a_{01}^2a_{20} + b_{10}^2b_{02} - 2b_{10}b_{01}b_{11} + b_{01}^2b_{20}, \\ & a_{10}^2b_{02} - 2a_{10}a_{01}b_{11} + a_{01}^2b_{20} + a_{20}b_{01}^2 - 2a_{11}b_{10}b_{01} + a_{02}b_{10}^2, \\ & \frac{a_{01}b_{01}}{a_{02} + b_{02}}, \quad \frac{a_{11}b_{02} - a_{02}b_{11}}{a_{02} - b_{02}}, \quad \frac{a_{10}b_{01} - a_{01}b_{10}}{a_{20}a_{02} - a_{11}^2 - b_{20}b_{02} + b_{11}^2}, \quad \frac{a_{02}^2 + b_{02}^2}{a_{02}b_{02}}, \\ & a_{10}a_{11}b_{01} - a_{10}a_{02}b_{10} - a_{10}b_{10}b_{02} + a_{10}b_{01}b_{11} - a_{01}a_{20}b_{01} + a_{01}a_{11}b_{10} + a_{01}b_{10}b_{11} - a_{01}b_{01}b_{20}. \end{aligned} $$

**Example 7.11 (Moment invariants for 2D images).** In this example, we will focus on moment invariants which is a popular tool in computer vision and pattern recognition [42, 43, 55, 59]. For our purposes, for any fixed integer $d$, the space of moments $M_d$ will simply be the space of nonconstant monomials of degree at most $d$ in variables $x$ and $y$. For example, $M_2$ is spanned by $x, y, x^2, xy, y^2$. The monomial $x^i y^j$ will be denoted by $m_{i,j}$. We will consider the action of the planar rotation group $SO(2)$ induced by its action on the $(x, y)$-plane. More precisely:

$$ \text{for } R = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \in \text{SO}(2) \text{ we define } R \cdot m(x, y) := m(ax + by, cx + dy). $$

<footnote>$^7$We decided to use this class of cross-sections based on the existing literature and private communication with Evelyne Hubert.</footnote>
For example, $m_{2,0} + m_{0,2} = x^2 + y^2$ is an invariant of this action, since rotations preserve the distance to the origin. The same is true for $m_{1,0}^2 + m_{0,1}^2$ which represents the same polynomial in $x$ and $y$ but is different as a function of the moments and, thus, is considered a different invariant.

The problem of computing invariants of this action goes back to the work of Hu [59] who proposed a set of seven invariants in the $d = 3$ case. This set turned out to be redundant [42] and was not complete (which is a desirable property, see [44]). A new irredundant set of invariants was proposed by Flusser [42], this set was not generating the whole field $k(M_3)^{SO(2)}$ but formed its transcendence basis. Another set was proposed in [55] which was still not a generating set of $k(M_3)^{SO(2)}$. In the same paper, the Hubert-Kogan algorithm was applied to the problem [55, p. 231] but the resulting invariants were deemed too complicated. We computed generators of $k(M_3)^{SO(2)}$ using the Hubert-Kogan algorithm with the cross-section $m_{1,1} = 0$ (the cross-section chosen using the same strategy as in Example 7.10). The generating set contained polynomials $m_{2,0} + m_{0,2}$ and $m_{2,0}m_{0,2} - m_{1,1}^2$, and 13 rational functions with degrees ranging from 8 to 10. After applying our simplification algorithm (with the result minimized using Algorithm 5) to this generating set, we arrive at the following set of seven polynomial generators:

$$f_1 = m_{2,0} + m_{0,2}, \quad f_2 = m_{2,0}m_{0,2} - m_{1,1}^2,$$

$$f_3 = m_{3,0}m_{1,2} - m_{2,1}^2 + m_{2,1}m_{0,3} - m_{1,2}^2, \quad f_4 = m_{3,0}^2 + 3m_{3,0}m_{1,2} + 3m_{2,1}m_{0,3} + m_{0,3}^2,$$

$$f_5 = m_{2,0}m_{2,1}m_{0,3} - m_{2,0}m_{1,2}^2 - m_{1,1}m_{3,0}m_{0,3} + m_{1,1}m_{2,1}m_{1,2} + m_{0,2}m_{3,0}m_{1,2} - m_{0,2}m_{2,1}^2,$$

$$f_6 = m_{2,0}m_{3,0}m_{0,3} - m_{2,0}m_{2,1}m_{1,2} - 2m_{1,1}m_{3,0}m_{1,2} + 2m_{1,1}m_{2,1}m_{0,3} - 2m_{1,1}m_{1,2}^2 \\ - m_{0,2}m_{3,0}m_{0,3} + m_{0,2}m_{2,1}m_{1,2},$$

$$f_7 = m_{2,0}m_{3,0}m_{1,2} + 2m_{2,0}m_{2,1}m_{0,3} + m_{2,0}m_{0,3}^2 - 2m_{1,1}m_{3,0}m_{2,1} - m_{1,1}m_{3,0}m_{0,3} - 3m_{1,1}m_{2,1}m_{1,2} \\ - 2m_{1,1}m_{1,2}m_{0,3} + m_{0,2}m_{3,0}^2 + 2m_{0,2}m_{3,0}m_{1,2} + m_{0,2}m_{2,1}m_{0,3}.$$

This generating set has lower degrees than the list from [42] (3 instead of 4) and smaller maximal monomial count than [55, p. 231] (10 instead of 15). On the other hand, our list recovers several classical invariants: $f_1$ and $f_2$ are equivalent to $\phi_1$ and $\phi_2$ in the Hu's list [59] following the notation from [42, p. 1405]; the same holds for $f_3$ and $f_4$ above and $\phi_3$ and $\phi_4$ ($\phi_3 = f_4 - 9f_3$ and $\phi_4 = f_4 - f_3$). This indicates that our automatic simplification procedure can produce simple generators coherent with the results in the pattern recognition literature. However, we would like to point out that, in this case, since $SO_2(\mathbb{C}) \cong \mathbb{C}^*$, more refined methods like [65] could be used instead for invariant computation (in the spirit of [42])$^8$.

## References

[1] S. Ahmed, N. Crepeau, P. R. Dessauer, A. Edozie, O. Garcia-Lopez, T. Grimsley, J. L. Garcia, V. Neri, and A. Shiu. Identifiability of directed-cycle and catenary linear compartmental models, 2026. URL http://dx.doi.org/10.1137/25M1728636.

[2] L. J. Allen. Some discrete-time SI, SIR, and SIS epidemic models. *Mathematical Biosciences*, 124 (1):83–105, Nov. 1994. URL http://dx.doi.org/10.1016/0025-5564(94)90025-6.

[3] E. S. Allman and J. A. Rhodes. *Mathematical Models in Biology: An Introduction*. Cambridge University Press, Oct. 2003. ISBN 9780511790911. URL http://dx.doi.org/10.1017/CBO9780511790911.

[4] A. Arnold and D. S. Roche. Multivariate sparse interpolation using randomized Kronecker substitutions. In *Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation*, ISSAC ’14, page 35–42, 2014. URL https://doi.org/10.1145/2608628.2608674.

[5] E. A. Arnold. Modular algorithms for computing Gröbner bases. *Journal of Symbolic Computation*, 35(4):403–419, 2003. URL https://doi.org/10.1016/S0747-7171(02)00140-2.

[6] G. Ballif, F. Clément, and R. Yvinec. Nonlinear compartmental modeling to monitor ovarian follicle population dynamics on the whole lifespan. *Journal of Mathematical Biology*, 89(9):43, 2024. URL https://doi.org/10.1007/s00285-024-02108-6.

[7] A. S. Bandeira, B. Blum-Smith, J. Kileel, J. Niles-Weed, A. Perry, and A. S. Wein. Estimation under group actions: Recovering orbits from invariants. *Applied and Computational Harmonic Analysis*, 66:236–319, Sept. 2023. URL http://dx.doi.org/10.1016/j.acha.2023.06.001.

<footnote n="8">⁸We thank Evelyne Hubert for this remark</footnote>
[8] X. R. Barreiro and A. F. Villaverde. On the origins and rarity of locally but not globally identifiable parameters in biological modeling. *IEEE Access*, 11:65457-65467, 2023. URL https://doi.org/10.1109/ACCESS.2023.3288998.

[9] A. Beauville. *The Lüroth Problem*, page 1–27. Springer International Publishing, 2016. ISBN 9783319462097. URL http://dx.doi.org/10.1007/978-3-319-46209-7_1.

[10] M. Ben-Or and P. Tiwari. A deterministic algorithm for sparse multivariate polynomial interpolation. In *Proceedings of the Twentieth Annual ACM Symposium on Theory of Computing*, STOC '88, page 301–309, 1988. URL https://doi.org/10.1145/62212.62241.

[11] P. D. Berk, R. B. Howe, J. R. Bloomer, and N. I. Berlin. Studies of bilirubin kinetics in normal adults. *Journal of Clinical Investigation*, 48(11):2176-2190, Nov. 1969. URL http://dx.doi.org/10.1172/JCI106184.

[12] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah. Julia: A fresh approach to numerical computing. *SIAM review*, 59(1):65–98, 2017. URL https://doi.org/10.1137/141000671.

[13] A. K. Binder. *Algorithms for Fields and an Application to a Problem in Computer Vision*. PhD thesis, Technische Universität München, 2009. URL https://mediatum.ub.tum.de/685465.

[14] F. Binder. Fast computations in the lattice of polynomial rational function fields. In *Proceedings of the 1996 international symposium on Symbolic and algebraic computation*, ISSAC '96, page 43–48. ACM Press, 1996. URL http://dx.doi.org/10.1145/236869.236895.

[15] B. Blum-Smith and H. Derksen. Generic orbits, normal bases, and generation degree for fields of rational invariants, 2025. URL https://arxiv.org/abs/2506.05650.

[16] D. K. Boku. *Gröbner Bases over Extention Fields of Q*. PhD thesis, Technische Universität Kaiserslautern, 2016. URL https://nbn-resolving.de/urn:nbn:de:hbz:386-kluedo-44289.

[17] D. K. Boku, W. Decker, and C. Fieker. `ffmodstd.lib`: A Singular version 4-0-3 library for computing Groebner bases of ideals in polynomial rings over algebraic function fields, 2016.

[18] H. Breloer. Rational invariants of even degree polynomials under the orthogonal group, 2025. URL https://arxiv.org/abs/2501.17504.

[19] M. Brickenstein. Slimgb: Gröbner Bases with Slim Polynomials. *Revista Matemática Complutense*, 23:453–466, 07 2010. URL https://doi.org/10.1007/s13163-009-0020-0.

[20] M. Bruno, J. Koschmieder, F. Wuest, P. Schaub, M. Fehling-Kaschek, J. Timmer, P. Beyer, and S. Al-Babili. Enzymatic study on AtCCD4 and AtCCD7 and their potential to form acyclic regulatory metabolites. *Journal of Experimental Botany*, 67(21):5993–6005, 10 2016. URL https://doi.org/10.1093/jxb/erw356.

[21] D. G. Cantor and H. Zassenhaus. A new algorithm for factoring polynomials over finite fields. *Mathematics of Computation*, 36:587–592, 1981. URL https://doi.org/10.2307/2007663.

[22] M. Capistran, M. Moreles, and B. Lara. Parameter estimation of some epidemic models. The case of recurrent epidemics caused by respiratory syncytial virus. *Bulletin of mathematical biology*, 71: 1890–901, 08 2009. URL https://doi.org/10.1007/s11538-009-9429-3.

[23] M. Castro, S. Ares, J. A. Cuesta, and S. Manrubia. The turning point and end of an expanding epidemic cannot be precisely forecast. *Proceedings of the National Academy of Sciences*, 117(42): 26190–26196, Oct. 2020. URL http://dx.doi.org/10.1073/pnas.2007868117.

[24] M. J. Chappell, K. R. Godfrey, and S. Vajda. Global identifiability of the parameters of nonlinear systems with specified inputs: A comparison of methods. *Mathematical Biosciences*, 102(1):41–73, 1990. URL https://doi.org/10.1016/0025-5564(90)90055-4.

[25] Q. Y. Chen, Z. Rapti, Y. Drossinos, J. Cuevas-Maraver, G. A. Kevrekidis, and P. G. Kevrekidis. Practical identifiability and parameter estimation of compartmental epidemiological models, 2024. URL https://arxiv.org/abs/2406.17827.

[26] C. Conradi and A. Shiu. Dynamics of posttranslational modification systems: Recent progress and future directions. *Biophysical Journal*, 114(3):507–515, 2018. URL https://doi.org/10.1016/j.bpj.2017.11.3787.
[27] F. Corniquel. Solving parametric polynomial systems using Generic Rational Univariate Representation. working paper or preprint, Jan. 2026. URL https://hal.science/hal-05487354.

[28] F. Crauste, J. Mafille, L. Boucinha, S. Djebali, O. Gandrillon, J. Marvel, and C. Arpin. Identification of nascent memory CD8 T cells and modeling of their ontogeny. *Cell Systems*, 4:306–317, 2017. ISSN 0022-5193. URL https://doi.org/10.1016/j.cels.2017.01.014.

[29] A. Cuyt and W.-s. Lee. Sparse interpolation of multivariate rational functions. *Theoretical Computer Science*, 412(16):1445–1456, 2011. URL https://doi.org/10.1016/j.tcs.2010.11.050.

[30] E. A. Dankwa, C. A. Donnelly, A. F. Brouwer, R. Zhao, M. P. Montgomery, M. K. Weng, and N. K. Martin. Estimating vaccination threshold and impact in the 2017–2019 hepatitis A virus outbreak among persons experiencing homelessness or who use drugs in Louisville, Kentucky, United States. *Vaccine*, 39(49):7182–7190, Dec. 2021. URL http://dx.doi.org/10.1016/j.vaccine.2021.10.001.

[31] E. A. Dankwa, A. F. Brouwer, and C. A. Donnelly. Structural identifiability of compartmental models for infectious disease transmission is influenced by data type. *Epidemics*, 41:100643, Dec. 2022. URL http://dx.doi.org/10.1016/j.epidem.2022.100643.

[32] S. Demignot and D. Domurado. Effect of prosthetic sugar groups on the pharmacokinetics of glucose-oxidase. *Drug Design and Delivery*, 1(4):333–348, May 1987. URL https://pubmed.ncbi.nlm.nih.gov/2855567/.

[33] R. A. Demillo and R. J. Lipton. A probabilistic remark on algebraic program testing. *Information Processing Letters*, 7(4):193–195, 1978. URL https://doi.org/10.1016/0020-0190(78)90067-4.

[34] A. Demin and S. Gowda. Groebner.jl: A package for Gröbner bases computations in Julia, 2023. URL https://arxiv.org/abs/2304.06935.

[35] A. Demin and J. van der Hoeven. Factoring sparse polynomials fast. *Journal of Complexity*, 88: 101934, 2025. URL https://doi.org/10.1016/j.jco.2025.101934.

[36] J. J. DiStefano III. *Dynamic Systems Biology Modeling and Simulation*. Academic Press Inc, 2015. ISBN 0124104118.

[37] R. Dong, C. Goodbrake, H. A. Harrington, and G. Pogudin. Differential elimination for dynamical models via projections with applications to structural identifiability. *SIAM Journal on Applied Algebra and Geometry*, 7(1):194–235, Mar. 2023. URL http://dx.doi.org/10.1137/22M1469067.

[38] R. Ehrenborg and G.-C. Rota. Apolarity and canonical forms for homogeneous polynomials. *European Journal of Combinatorics*, 14(3):157–181, 1993. URL https://doi.org/10.1006/eujc.1993.1022.

[39] J.-C. Faugére. A new efficient algorithm for computing Gröbner bases (F4). *Journal of Pure and Applied Algebra*, 139(1):61–88, 1999. URL https://doi.org/10.1016/S0022-4049(99)00005-5.

[40] C. Fieker, W. Hart, T. Hofmann, and F. Johansson. Nemo/Hecke: Computer algebra and number theory packages for the Julia programming language. In *Proceedings of the 2017 ACM on International Symposium on Symbolic and Algebraic Computation*, ISSAC '17, pages 157–164, 2017. URL https://doi.acm.org/10.1145/3087604.3087611.

[41] L. Fink. Model selection in deterministic models of mRNA transfection, Master's thesis, 2015.

[42] J. Flusser. On the independence of rotation moment invariants. *Pattern Recognition*, 33 (9):1405–1410, Sept. 2000. URL http://dx.doi.org/10.1016/S0031-3203(99)00127-2.

[43] J. Flusser, T. Suk, and B. Zitová. *Moments and Moment Invariants in Pattern Recognition*. Wiley, Oct. 2009. ISBN 9780470684757. URL http://dx.doi.org/10.1002/9780470684757.

[44] J. Flusser, T. Suk, and B. Zitová. Complete and incomplete sets of invariants. *Journal of Mathematical Imaging and Vision*, 63(7):917–922, May 2021. URL http://dx.doi.org/10.1007/s10851-021-01039-x.
[45] A. Fokas, J. Cuevas-Maraver, and P. Kevrekidis. A quantitative framework for exploring exit strategies from the COVID-19 lockdown. *Chaos, Solitons & Fractals*, 140:110244, 2020. URL http://dx.doi.org/10.1016/j.chaos.2020.110244.

[46] K. A. Fujita, Y. Toyoshima, S. Uda, Y. I. Ozaki, H. Kubota, and S. Kuroda. Decoupling of receptor and downstream signals in the Akt pathway by its low-pass filter characteristics. *Science Signaling*, 3(132):ra56-ra56, 2010. URL https://doi.org/10.1126/scisignal.2000810.

[47] J. L. Gevertz, J. M. Greene, C. H. Sanchez-Tapia, and E. D. Sontag. A novel COVID-19 epidemiological model with explicit susceptible and asymptomatic isolation compartments reveals unexpected consequences of timing social distancing. *Journal of Theoretical Biology*, 510:110539, 2021. URL https://doi.org/10.1016/j.jtbi.2020.110539.

[48] M. Giesbrecht and D. S. Roche. Diversification improves interpolation. In *Proceedings of the 36th International Symposium on Symbolic and Algebraic Computation*, ISSAC '11, page 123–130, 2011. ISBN 9781450306751. URL https://doi.org/10.1145/1993886.1993909.

[49] B. C. Goodwin. Oscillatory behavior in enzymatic control processes. *Advances in Enzyme Regulation*, 3:425–437, 1965. URL https://doi.org/10.1016/0065-2571(65)90067-1.

[50] P. Görlach, E. Hubert, and T. Papadopoulo. Rational invariants of even ternary forms under the orthogonal group. *Foundations of Computational Mathematics*, 19:1315–1361, 2019. URL https://doi.org/10.1007/s10208-018-9404-1.

[51] H.-G. Gräbe. The SymbolicData Benchmark Problems Collection of Polynomial Systems. In *Proceedings of the Workshop on Under- and Overdetermined Systems of Algebraic or Differential Equations*, pages 57–76. IAS, Univ. Karlsruhe, 2002. URL https://symbolicdata.github.io/Papers/karlsruhe-02.pdf.

[52] T. R. Grandjean, M. J. Chappell, J. W. Yates, and N. D. Evans. Structural identifiability analyses of candidate models for in vitro Pitavastatin hepatic uptake. *Computer Methods and Programs in Biomedicine*, 114(3):e60–e69, 2014. URL https://doi.org/10.1016/j.cmpb.2013.06.013.

[53] B. Grenet, J. van der Hoeven, and G. Lecerf. Deterministic root finding over finite fields using Graeffe transforms. *Applicable Algebra in Engineering, Communication and Computing*, 27:237–257, 06 2016. URL https://doi.org/10.1007/s00200-015-0280-5.

[54] J. Gutierrez, R. Rubio, and D. Sevilla. Unirational fields of transcendence degree one and functional decomposition. In *Proceedings of the 2001 international symposium on Symbolic and algebraic computation*, ISSAC01, page 167–174. ACM, July 2001. URL http://dx.doi.org/10.1145/384101.384124.

[55] M. S. Hickman. Geometric moments and their invariants. *Journal of Mathematical Imaging and Vision*, 44(3):223–235, Dec. 2011. URL http://dx.doi.org/10.1007/s10851-011-0323-x.

[56] H. Hong, A. Ovchinnikov, G. Pogudin, and C. Yap. SIAN: software for structural identifiability analysis of ODE models. *Bioinformatics*, 35(16):2873–2874, 2019. URL http://dx.doi.org/10.1093/bioinformatics/bty1069.

[57] H. Hong, A. Ovchinnikov, G. Pogudin, and C. Yap. Global identifiability of differential models. *Communications on Pure and Applied Mathematics*, 73(9):1831–1879, 2020. URL https://doi.org/10.1002/cpa.21921.

[58] J. Hu and M. Monagan. A fast parallel sparse polynomial gcd algorithm. In *Proceedings of the 2016 ACM International Symposium on Symbolic and Algebraic Computation*, ISSAC '16, page 271–278, 2016. ISBN 9781450343800. URL https://doi.org/10.1145/2930889.2930903.

[59] M.-K. Hu. Visual pattern recognition by moment invariants. *IEEE Transactions on Information Theory*, 8(2):179–187, Feb. 1962. ISSN 0018-9448. URL http://dx.doi.org/10.1109/TIT.1962.1057692.

[60] E. Hubert. The AIDA Maple package. URL https://www-sop.inria.fr/members/Evelyne.Hubert/aida/. accessed on November 3, 2025.
[61] E. Hubert and M. Jalard. Rationality of the invariant field for a class of representations of the real orthogonal groups. *Journal of Algebra*, 682:109–130, 2025. URL http://dx.doi.org/10.1016/j.jalgebra.2025.05.033.

[62] E. Hubert and M. Jalard. Orbit separation and stratification by isotropy classes of piezoelectricity tensors. *Journal of Pure and Applied Algebra*, 229(9):108034, 2025. URL https://doi.org/10.1016/j.jpaa.2025.108034.

[63] E. Hubert and M. Jalard. Algebraically independent generators for the invariant field of $so_3$ and $o_3$ representations $r^3 \oplus h$. to appear in *SIAM Journal on Applied Algebra and Geometry*, (1), 2026. URL https://inria.hal.science/hal-04604969.

[64] E. Hubert and I. A. Kogan. Rational invariants of a group action. Construction and rewriting. *Journal of Symbolic Computation*, 42(1-2):203–217, Jan. 2007. URL http://dx.doi.org/10.1016/j.jsc.2006.03.005.

[65] E. Hubert and G. Labahn. Scaling invariants and symmetry reduction of dynamical systems. *Foundations of Computational Mathematics*, 13(4):479–516, Aug. 2013. URL http://dx.doi.org/10.1007/s10208-013-9165-9.

[66] N. Idrees, G. Pfister, and S. Steidel. Parallelization of modular algorithms. *Journal of Symbolic Computation*, 46(6):672–684, 2011. URL https://doi.org/10.1016/j.jsc.2011.01.003.

[67] I. Ilmer, A. Ovchinnikov, and G. Pogudin. *Web-Based Structural Identifiability Analyzer*, page 254–265. Springer International Publishing, 2021. ISBN 9783030856335. URL http://dx.doi.org/10.1007/978-3-030-85633-5_17.

[68] I. Ilmer, A. Ovchinnikov, G. Pogudin, and P. Soto. More efficient identifiability verification in ODE models by reducing non-identifiability, 2022. URL https://arxiv.org/abs/2204.01623.

[69] A. Jiménez-Pastor and G. Pogudin. Computing exact nonlinear reductions of dynamical models. *ACM Communications in Computer Algebra*, 56(2):25–31, June 2022. URL http://dx.doi.org/10.1145/3572867.3572869.

[70] A. Jinadu and M. Monagan. Solving parametric linear systems using sparse rational function interpolation. In *Computer Algebra in Scientific Computing*, CASC 2023, page 233–254, Berlin, Heidelberg, 2023. Springer-Verlag. ISBN 978-3-031-41723-8. URL https://doi.org/10.1007/978-3-031-41724-5_13.

[71] T. A. Joseph, L. Shenhav, J. B. Xavier, E. Halperin, and I. Pe'er. Compositional Lotka-Volterra describes microbial dynamics in the simplex. *PLOS Computational Biology*, 16(5):1–22, 05 2020. URL https://doi.org/10.1371/journal.pcbi.1007917.

[72] A. Joux and V. Vitse. A variant of the F4 algorithm. In A. Kiayias, editor, *Topics in Cryptology – CT-RSA 2011*, pages 356–375, Berlin, Heidelberg, 2011. Springer Berlin Heidelberg. ISBN 978-3-642-19074-2. URL https://doi.org/10.1007/978-3-642-19074-2_23.

[73] E. Kaltofen and W.-s. Lee. Early termination in sparse interpolation algorithms. *Journal of Symbolic Computation*, 36(3):365–400, 2003. URL https://doi.org/10.1016/S0747-7171(03)00088-9.

[74] E. Kaltofen and L. Yagati. Improved sparse multivariate polynomial interpolation algorithms. In P. Gianni, editor, *Symbolic and Algebraic Computation*, pages 467–474, Berlin, Heidelberg, 1989. Springer Berlin Heidelberg. URL https://doi.org/10.1007/3-540-51084-2_44.

[75] E. Kaltofen, W.-s. Lee, and A. A. Lobo. Early termination in Ben-Or/Tiwari sparse interpolation and a hybrid of Zippel's algorithm. In *Proceedings of the 2000 International Symposium on Symbolic and Algebraic Computation*, page 192–201, 2000. URL https://doi.org/10.1145/345542.345629.

[76] Y. Kawano and T. Ohtsuka. An algebraic approach to local observability at an initial state for discrete-time polynomial systems. *IFAC Proceedings Volumes*, 44(1):6449–6453, Jan. 2011. URL http://dx.doi.org/10.3182/20110828-6-IT-1002.00336.

[77] G. Kemper. A constructive approach to Noether's problem. *Manuscripta Mathematica*, 90(1): 343–363, Dec. 1996. URL http://dx.doi.org/10.1007/BF02568311.
[78] G. Kemper. The computation of invariant fields and a constructive version of a theorem by Rosenlicht. *Transformation Groups*, 12(4):657–670, Nov. 2007. URL http://dx.doi.org/10.1007/s00031-007-0056-5.

[79] S. Khodadad and M. Monagan. Fast rational function reconstruction. In *Proceedings of the 2006 International Symposium on Symbolic and Algebraic Computation*, ISSAC '06, page 184–190, 2006. URL https://doi.org/10.1145/1145768.1145801.

[80] M. Kim and F. Milner. A mathematical model of epidemics with screening and variable infectivity. *Mathematical and Computer Modelling*, 21(7):29–42, 1995. URL https://doi.org/10.1016/0895-7177(95)00029-2.

[81] I. A. Kogan. Invariants: Computation and applications. In *Proceedings of the 2023 International Symposium on Symbolic and Algebraic Computation*, ISSAC 2023, page 31–40. ACM, 2023. URL http://dx.doi.org/10.1145/3597066.3597149.

[82] A. Kumar and P. Daoutidis. *Control of nonlinear differential algebraic equation systems with applications to chemical processes*. Chapman & Hall/CRC Research Notes in Mathematics Series. CRC Press, Boca Raton, FL, Feb. 1999.

[83] S. Kuroda. A counterexample to the fourteenth problem of Hilbert in dimension three. *Michigan Mathematical Journal*, 53(1), 2005. URL http://dx.doi.org/10.1307/MMJ/1114021089.

[84] S. Lang. *Algebra*. Springer, 3rd edition, 2002. URL https://doi.org/10.1007/978-1-4613-0041-0.

[85] Y. Lecourtier, F. Lamnabhi-Lagarrigue, and E. Walter. A method to prove that nonlinear models can be unidentifiable. In *26th IEEE Conference on Decision and Control*, volume 26, pages 2144–2145, 1987. URL https://doi.org/10.1109/CDC.1987.272467.

[86] E. C. Lee, M. R. Kelly, B. M. Ochocki, S. M. Akinwumi, K. E. Hamre, J. H. Tien, and M. C. Eisenberg. Model distinguishability and inference robustness in mechanisms of cholera transmission and loss of immunity. *Journal of Theoretical Biology*, 420:68–81, 2017. URL https://doi.org/10.1016/j.jtbi.2017.01.032.

[87] T. S. Ligon, F. Fröhlich, O. T. Chiş, J. R. Banga, E. Balsa-Canto, and J. Hasenauer. GenSSI 2.0: multi-experiment structural identifiability analysis of SBML models. *Bioinformatics*, 34(8):1421–1423, 11 2017. URL https://doi.org/10.1093/bioinformatics/btx735.

[88] C. Lyzell, T. Glad, M. Enqvist, and L. Ljung. Difference algebra and system identification. *Automatica*, 47(9):1896–1904, 2011. URL https://doi.org/10.1016/j.automatica.2011.06.013.

[89] A. K. Manrai and J. Gunawardena. The geometry of multisite phosphorylation. *Biophysical Journal*, 95(12):5533–5543, 2008. URL https://doi.org/10.1529/biophysj.108.140632.

[90] G. Massonis, J. R. Banga, and A. F. Villaverde. Structural identifiability and observability of compartmental models of the COVID-19 pandemic. *Annual Reviews in Control*, 51:441–459, 2021. URL https://doi.org/10.1016/j.arcontrol.2020.12.001.

[91] N. Meshkat, Z. Rosen, and S. Sullivant. Algebraic tools for the analysis of state space models. In *The 50th Anniversary of Gröbner Bases*, page 171–205. Mathematical Society of Japan. URL http://dx.doi.org/10.2969/aspm/07710171.

[92] N. Meshkat, M. Eisenberg, and J. J. DiStefano III. An algorithm for finding globally identifiable parameter combinations of nonlinear ode models using Gröbner bases. *Mathematical Biosciences*, 222(2):61–72, 2009. URL https://doi.org/10.1016/j.mbs.2009.08.010.

[93] N. Meshkat, C. Anderson, and J. J. DiStefano III. Finding identifiable parameter combinations in nonlinear ODE models and the rational reparameterization of their input–output equations. *Mathematical Biosciences*, 233(1):19–31, 2011. URL http://dx.doi.org/10.1016/j.mbs.2011.06.001.

[94] N. Meshkat, C. Kuo, and J. J. DiStefano III. On finding and using identifiable parameter combinations in nonlinear dynamic systems biology models and COMBOS: A novel web implementation. *PloS one*, 9:e110261, 10 2014. URL https://doi.org/10.1371/journal.pone.0110261.
[95] J. S. Milne. *Fields and Galois Theory*. Kea Books, Ann Arbor, MI, 2022.

[96] P. Moate, R. Boston, T. Jenkins, and I. Lean. Kinetics of ruminal lipolysis of triacylglycerol and biohydrogenation of long-chain fatty acids: New insights from old data. *Journal of Dairy Science*, 91(2):731–742, 2008. URL https://doi.org/10.3168/jds.2007-0398.

[97] M. Monagan and R. Pearce. An algorithm for splitting polynomial systems based on F4. In *Proceedings of the International Workshop on Parallel Symbolic Computation*, PASCO 2017, 2017. URL https://doi.org/10.1145/3115936.3115948.

[98] Y. Mukhina and G. Pogudin. Projecting dynamical systems via a support bound, 2025. URL https://arxiv.org/abs/2501.13680.

[99] J. Müller-Quade and T. Beth. Calculating generators for invariant fields of linear algebraic groups. In *Applied algebra, algebraic algorithms and error-correcting codes* (Honolulu, HI, 1999), volume 1719 of *Lecture Notes in Computer Science*. Springer, 1999. URL https://doi.org/10.1007/3-540-46796-3_37.

[100] J. Müller-Quade and R. Steinwandt. Gröbner bases applied to finitely generated field extensions. *Journal of Symbolic Computation*, 30(4):469–490, Oct. 2000. URL http://dx.doi.org/10.1006/jsco.1999.0417.

[101] J. Müller-Quade and R. Steinwandt. Basic algorithms for rational function fields. *Journal of Symbolic Computation*, 27(2):143–170, 1999. URL https://doi.org/10.1006/jsco.1998.0246.

[102] M. Nagata. On the 14-th problem of Hilbert. *American Journal of Mathematics*, 81(3):766, 1959. URL http://dx.doi.org/10.2307/2372927.

[103] C. Nakul. Introduction to SEIR models, 2017. URL https://indico.ictp.it/event/7960/session/3/contribution/19/material/slides/0.pdf.

[104] F. Ollivier. *Le problème de l'identifiabilité structurelle globale: approche théorique, méthodes effectives et bornes de complexité*. PhD thesis, École polytechnique, 1990. URL https://www.theses.fr/1990EPXX0009.

[105] A. Ovchinnikov, A. Pillay, G. Pogudin, and T. Scanlon. Computing all identifiable functions of parameters for ODE models. *Systems & Control Letters*, 157:105030, 2021. URL https://doi.org/10.1016/j.sysconle.2021.105030.

[106] A. Ovchinnikov, G. Pogudin, and P. Thompson. Parameter identifiability and input-output equations. *Applicable Algebra in Engineering, Communication and Computing*, 34(2):165–182, 2021. URL https://doi.org/10.1007/s00200-021-00486-8.

[107] G. Pogudin. *Symbolic Transformations of Dynamical Models*. Habilitation à diriger des recherches, Institute Polytechnique de Paris, 2024. URL https://hal.science/tel-05041736.

[108] R. Prony. Essai experimental et analytique sur les lois de la dilatabilite de fluides elastiques et sur celles da la force expansion de la vapeur de l'alcool, a differentes temperatures. *J. de l'École Polytechnique Floréal et Plairial, an III*, 22:24–76, 1795.

[109] V. Raia, M. Schilling, M. Böhm, B. Hahn, A. Kowarsch, A. Raue, C. Sticht, S. Bohl, M. Saile, P. Möller, N. Gretz, J. Timmer, F. Theis, W.-D. Lehmann, P. Lichter, and U. Klingmüller. Dynamic mathematical modeling of IL13-induced signaling in hodgkin and primary mediastinal B-cell lymphoma allows prediction of therapeutic targets. *Cancer Research*, 71(3):693–704, 2011. URL https://doi.org/10.1158/0008-5472.CAN-10-2987.

[110] C. H. Remien, M. J. Eckwright, and B. J. Ridenhour. Structural identifiability of the generalized Lotka–Volterra model for microbiome studies. *Royal Society Open Science*, 8(7):20378, 2021. URL https://doi.org/10.1098/rsos.201378.

[111] X. Rey Barreiro and A. F. Villaverde. Benchmarking tools for a priori identifiability analysis. *Bioinformatics*, 39(2), Jan. 2023. URL http://dx.doi.org/10.1093/bioinformatics/btad065.

[112] L. Robbiano and M. Sweedler. *Subalgebra bases*, page 61–87. Springer Berlin Heidelberg, 1990. ISBN 9783540471363. URL http://dx.doi.org/10.1007/BFb0085537.
[113] K. Roosa and G. Chowell. Assessing parameter identifiability in compartmental dynamic models using a computational approach: application to infectious disease transmission models. *Theoretical Biology and Medical Modelling*, 16, 01 2019. URL https://doi.org/10.1186/s12976-018-0097-6.

[114] M. P. Saccomani and G. Bellu. DAISY: An efficient tool to test global identifiability. some case studies. In *2008 16th Mediterranean Conference on Control and Automation*, pages 1723–1728, 2008. URL https://doi.org/10.1109/MED.2008.4602152.

[115] M. P. Saccomani, S. Audoly, G. Bellu, and L. D'Angiò. Examples of testing global identifiability of biological and biomedical models with the DAISY software. *Computers in Biology and Medicine*, 40 (4):402–407, 2010. URL https://doi.org/10.1016/j.compbiomed.2010.02.004.

[116] T. Sauer, T. Berry, D. Ebeigbe, M. M. Norton, A. J. Whalen, and S. J. Schiff. Identifiability of infection model parameters early in an epidemic. *SIAM Journal on Control and Optimization*, 60(2): S27–S48, Nov. 2021. URL http://dx.doi.org/10.1137/20M1353289.

[117] J. T. Schwartz. Fast probabilistic algorithms for verification of polynomial identities. *J. ACM*, 27(4): 701–717, Oct. 1980. URL https://doi.org/10.1145/322217.322225.

[118] A. Sedoglavic. A probabilistic algorithm to test local algebraic observability in polynomial time. *Journal of Symbolic Computation*, 33(5):735–755, 2002. URL https://doi.org/10.1006/jsco.2002.0532.

[119] C. Sin. Algorithmic parameter space reduction of a systems biology model: A case study, 2012. URL https://escholarship.org/uc/item/4r25f875. Master thesis.

[120] K. E. Smith and J. Rosenberg. Rational and non-rational algebraic varieties: Lectures of János Kollár, 1997. URL https://arxiv.org/abs/alg-geom/9707013.

[121] M. Sweedler. *Using Gröbner bases to determine the algebraic and transcendental nature of field extensions: Return of the killer tag variables*, page 66–75. Springer Berlin Heidelberg, 1993. ISBN 9783540476306. URL http://dx.doi.org/10.1007/3-540-56686-4_34.

[122] C. Traverso. Gröbner trace algorithms. In P. Gianni, editor, *Symbolic and Algebraic Computation*, pages 125–138, Berlin, Heidelberg, 1989. Springer Berlin Heidelberg. URL https://doi.org/10.1007/3-540-51084-2_12.

[123] N. Tuncer and T. T. Le. Structural and practical identifiability analysis of outbreak models. *Mathematical Biosciences*, 299:1–18, 2018. URL http://dx.doi.org/10.1016/j.mbs.2018.02.004.

[124] J. van der Hoeven and G. Lecerf. On sparse interpolation of rational functions and gcds. *ACM Commun. Comput. Algebra*, 55(1):1–12, May 2021. URL https://doi.org/10.1145/3466895.3466896.

[125] J. van der Hoeven and G. Lecerf. Sparse polynomial interpolation: faster strategies over finite fields. *Applicable Algebra in Engineering, Communication and Computing*, apr 2024. URL https://doi.org/10.1007/s00200-024-00655-5.

[126] J. von zur Garthen and J. Gerhard. *Modern Computer Algebra*. Cambridge University Press, 3rd edition, 2013.

[127] P. S. Wang, M. J. T. Guy, and J. H. Davenport. P-adic reconstruction of rational numbers. *SIGSAM Bull.*, 16(2):2–3, May 1982. URL https://doi.org/10.1145/1089292.1089293.

[128] H. Weyl. *The Classical Groups: Their Invariants and Representations*. Princeton University Press, 2nd edition, 1997.

[129] B. Whipple, T. A. Miura, and E. A. Hernandez-Vargas. Modeling the CD8+ T cell immune response to influenza infection in adult and aged mice. *Journal of Theoretical Biology*, 593:111898, 2024. URL https://doi.org/10.1016/j.jtbi.2024.111898.

[130] D. Wodarz and M. A. Nowak. Specific therapy regimes could lead to long-term immunological control of HIV. *Proceedings of the National Academy of Sciences*, 96(25):14464–14469, 1999. URL https://www.pnas.org/doi/abs/10.1073/pnas.96.25.14464.
[131] W.-t. Zha, F.-r. Pang, N. Zhou, B. Wu, Y. Liu, Y.-b. Du, X.-q. Hong, and Y. Lv. Research about the optimal strategies for prevention and control of varicella outbreak in a school in a central city of China: based on an SEIR dynamic model. *Epidemiology and Infection*, 148:e56, 2020. URL https://doi.org/10.1017/S0950268819002188.

[132] W. Zheng. Total variation regularization for compartmental epidemic models with time-varying dynamics, 2020. URL https://arxiv.org/abs/2004.00412.

[133] R. Zippel. Probabilistic algorithms for sparse polynomials. In *Proceedings of the International Symposium on Symbolic and Algebraic Computation*, EUROSAM '79, page 216–226, Berlin, Heidelberg, 1979. Springer-Verlag. ISBN 3540095195.
# A Appendix

This section lists our benchmark examples that arise from the problem of structural identifiability of ODE models<fnref n="9" /> (see Section 7.1). Many of these examples have been collected in [8, 37, 56, 87, 118]. For each example, we report:

* information about the original generating set;

* the names of the indeterminates;

* the original generating set (when it is not too large);

* the result of our algorithm, that is, a simplified generating set.

We also provide the generating sets in plain text format, together with the code to produce them<fnref n="10" />. The files that exceed 1 MB in size are not provided, but code in Julia is provided to reproduce the files.

## 1. Bilirubin [94, Eq. (14)].

Original generating set information: 7 indeterminates; 7 non-constant functions; maximal total degrees of numerator and denominator are (4, 0); 390 bytes total in string representation.

Indeterminates: $k_{01}, k_{12}, k_{13}, k_{14}, k_{21}, k_{31}, k_{41}$.

Original generating set:

$$ 
\begin{aligned}
-k_{12}k_{13}k_{14}, & k_{01}k_{12}k_{13}k_{14}, -k_{12} - k_{13} - k_{14}, \\
-k_{12}k_{13} - k_{12}k_{14} - k_{13}k_{14}, & k_{01} + k_{12} + k_{13} + k_{14} + k_{21} + k_{31} + k_{41}, \\
k_{01}k_{12}k_{13} + k_{01}k_{12}k_{14} + k_{01}k_{13}k_{14} + k_{12}k_{13}k_{14} + k_{12}k_{13}k_{41} + k_{12}k_{14}k_{31} + k_{13}k_{14}k_{21}, \\
k_{01}k_{12} + k_{01}k_{13} + k_{01}k_{14} + k_{12}k_{13} + k_{12}k_{14} + k_{12}k_{31} + k_{12}k_{41} + k_{13}k_{14} + k_{13}k_{21} + k_{13}k_{41} + k_{14}k_{21} + k_{14}k_{31}.
\end{aligned}
$$

Result of our algorithm:

$$ 
\begin{aligned}
k_{01}, k_{12}k_{13}k_{14}, k_{21}k_{31}k_{41}, k_{12} + k_{13} + k_{14}, k_{21} + k_{31} + k_{41}, & k_{12}k_{13} + k_{12}k_{14} + k_{13}k_{14}, \\
k_{21}k_{31} + k_{21}k_{41} + k_{31}k_{41}, & k_{12}k_{31} + k_{12}k_{41} + k_{13}k_{21} + k_{13}k_{41} + k_{14}k_{21} + k_{14}k_{31}.
\end{aligned}
$$

Result is algebraically independent over $\mathbb{C}$: no.

## 2. Biohydrogenation [96, Eq. (6)-(9)].

Original generating set information: 6 indeterminates; 37 non-constant functions; maximal total degrees of numerator and denominator are (6, 1); 1.2 KB total in string representation.

Indeterminates: $k_5, k_6, k_7, k_8, k_9, k_{10}$.

Original generating set: too large to be listed.

Result of our algorithm:

$$ k_5, k_6, k_7, k_9^2, k_9k_{10}, 2k_8 + k_{10}. $$

Result is algebraically independent over $\mathbb{C}$: yes.

## 3. Bruno2016 [20].

Original generating set information: 6 indeterminates; 4 non-constant functions; maximal total degrees of numerator and denominator are (2, 0); 166 bytes total in string representation.

Indeterminates: $k_{\text{zea}}, k_{\beta}, k_{\text{cry-OH}}, k_{\beta-10}, k_{\text{cry-\beta}}, k_{\text{OH-\beta-10}}$.

Original generating set:

$$ k_{\beta}, k_{\beta-10} + k_{\text{cry-OH}} + k_{\text{cry-\beta}}, k_{\beta-10}k_{\text{cry-OH}} + k_{\beta-10}k_{\text{cry-\beta}}, k_{\beta}^2 - k_{\beta}k_{\text{cry-OH}} - k_{\beta}k_{\text{cry-\beta}}. $$

Result of our algorithm:

$$ k_{\beta}, k_{\beta-10}, k_{\text{cry-OH}} + k_{\text{cry-\beta}}. $$

Result is algebraically independent over $\mathbb{C}$: yes.

## 4. CGV1990 [24, 118].

Original generating set information: 9 indeterminates; 412 non-constant functions; maximal total degrees of numerator and denominator are (18, 6); 90.0 KB total in string representation.

Indeterminates: $R, S, V_3, k_3, k_4, k_5, k_6, k_7, V_{36}$.

<footnote n="9">⁹The sources of the ODE models are available at https://github.com/SciML/StructuralIdentifiability.jl/blob/v0.5.18/benchmarking/benchmarks.jl</footnote>

<footnote n="10">¹⁰The data in plain text format is available at https://github.com/pogudingleb/RationalFunctionFields.jl/tree/f7d4943aa7432bd10c206b59f63909aca52c909d/paper/data</footnote>
Original generating set: too large to be listed.

Result of our algorithm:

$$k_3, k_4, k_6, k_7, RV_3 + SV_{36}, RV_{36} + \frac{1}{25}SV_3, \frac{V_3k_5 + 5V_{36}k_5}{V_3}, RV_{36}k_5 + \frac{1}{5}SV_{36}k_5.$$

Result is algebraically independent over $\mathbb{C}$: yes.

5. CRN [26, Eq. (2)].

Original generating set information: 6 indeterminates; 69 non-constant functions; maximal total degrees of numerator and denominator are (5, 2); 2.0 KB total in string representation.

Indeterminates: $k_1, k_2, k_3, k_4, k_5, k_6$.

Original generating set: too large to be listed.

Result of our algorithm:

$$k_1, k_2, k_3, k_4, k_5, k_6.$$

Result is algebraically independent over $\mathbb{C}$: yes.

6. Covid1 [47, Eq. (8)-(13)].

Original generating set information: 10 indeterminates; 302 non-constant functions; maximal total degrees of numerator and denominator are (16, 4); 214.5 KB total in string representation.

Indeterminates: $f, \beta_A, \beta_I, \epsilon_A, \epsilon_S, h_1, h_2, \gamma_{AI}, \gamma_{IR}, \delta$.

Original generating set: too large to be listed.

Result of our algorithm:

$$f, \beta_A, \beta_I, \epsilon_A, \epsilon_S, h_1, h_2, \gamma_{AI}, \delta + \gamma_{IR}.$$

Result is algebraically independent over $\mathbb{C}$: yes.

7. Covid3 [103, Slide 35].

Original generating set information: 8 indeterminates; 7 non-constant functions; maximal total degrees of numerator and denominator are (4, 2); 247 bytes total in string representation.

Indeterminates: $k, \Lambda, N, b, e, g, m, r$.

Original generating set:

$$\frac{br}{kN}, e + g + 3m, -e - g - 2m, em + gm + 2m^2, \\ \frac{ber + bgr + 2bmr}{kN}, \frac{begr + bemr + bgmr + bm^2r}{kN}, \frac{-\Lambda be + Negm + Nem^2 + Ngm^2 + Nm^3}{N}.$$

Result of our algorithm:

$$m, eg, e + g, \frac{k\Lambda}{gr}, \frac{kN}{br}.$$

Result is algebraically independent over $\mathbb{C}$: yes.

8. Crauste [28, Figure 3C].

Original generating set information: 13 indeterminates; 30 non-constant functions; maximal total degrees of numerator and denominator are (4, 2); 1.2 KB total in string representation.

Indeterminates: $\mu_M, \mu_N, \mu_P, \mu_{EE}, \mu_{LE}, \mu_{LL}, \mu_{PE}, \mu_{PL}, \rho_E, \rho_P, \delta_{EL}, \delta_{LM}, \delta_{NE}$.

Original generating set: too large to be listed.

Result of our algorithm:

$$\mu_M, \mu_N, \mu_P, \mu_{EE}, \mu_{LE}, \mu_{LL}, \mu_{PE}, \mu_{PL}, \delta_{EL}, \delta_{LM}, \frac{\rho_E}{\rho_P}, \frac{\delta_{NE}}{\rho_P}.$$

Result is algebraically independent over $\mathbb{C}$: yes.

9. EAIHRD [45, Eq. (1)-(6)], equations copied from [25, Eq. (6a)-(6f)].

Original generating set information: 10 indeterminates; 316 non-constant functions; maximal total degrees of numerator and denominator are (22, 6); 11.0 MB total in string representation.

Indeterminates: $N, a, d, h, s, c_1, c_2, r_1, r_2, r_3$.

Original generating set: too large to be listed.

Result of our algorithm:
$$r_1, d+r_3, a+h+r_2+s, ah+ar_2+hs+r_2s, \frac{dhs}{ac_1+c_2s}, \frac{ac_1h+ac_1r_2+c_2r_1s}{ac_1+c_2s}.$$

Result is algebraically independent over $\mathbb{C}$: yes.

10. Fujita [46], [37, Example Akt pathway].

Original generating set information: 16 indeterminates; 217 non-constant functions; maximal total degrees of numerator and denominator are (9, 5); 24.9 KB total in string representation.

Indeterminates: $a_1, a_2, a_3, r_{11}, r_{12}, r_{21}, r_{22}, r_{31}, r_{41}, r_{51}, r_{52}, r_{61}, r_{71}, r_{81}, r_{91}$.

Original generating set: too large to be listed.

Result of our algorithm:

$$r_{22}, r_{31}, r_{41}, r_{52}, r_{61}, r_{71}, r_{81}, r_{11} - r_{12} - r_{91}, \frac{a_1}{r_{51}}, \frac{a_2}{r_{51}}, \frac{a_3}{r_{51}}, \frac{r_{21}}{r_{51}}.$$

Result is algebraically independent over $\mathbb{C}$: yes.

11. Goodwin [49], [56, Example Goodwin oscillator in Supplementary materials].

Original generating set information: 7 indeterminates; 89 non-constant functions; maximal total degrees of numerator and denominator are (12, 0); 4.6 KB total in string representation.

Indeterminates: $b, c, \beta, \gamma, \alpha, \delta, \sigma$.

Original generating set: too large to be listed.

Result of our algorithm:

$$b, c, \sigma, \beta\delta, \beta + \delta.$$

Result is algebraically independent over $\mathbb{C}$: yes.

12. HighDimNonLin [115, Section 6].

Original generating set information: 22 indeterminates; 42 non-constant functions; maximal total degrees of numerator and denominator are (2, 0); 271 bytes total in string representation.

Indeterminates: $k_m, p_1, p_2, p_3, p_4, p_5, p_6, p_7, p_8, p_9, v_m, p_{10}, p_{11}, p_{12}, p_{13}, p_{14}, p_{15}, p_{16}, p_{17}, p_{18}, p_{19}, p_{20}$.

Original generating set:

$$k_m, p_1, p_2, p_3, p_4, p_5, p_6, p_7, p_8, p_9, -k_m, -p_1, -p_2, -p_3, -p_4, -p_5, -p_6, -p_7, -p_8, -p_9, p_{10}, p_{11}, p_{12}, p_{13}, p_{14}, p_{15}, p_{16}, p_{17}, p_{18}, p_{19}, p_{20}, -p_{10}, -p_{11}, -p_{12}, -p_{13}, -p_{14}, -p_{15}, -p_{16}, -p_{17}, -p_{18}, -p_{19}, k_m p_1 + v_m.$$

Result of our algorithm:

$$k_m, p_1, p_2, p_3, p_4, p_5, p_6, p_7, p_8, p_9, v_m, p_{10}, p_{11}, p_{12}, p_{13}, p_{14}, p_{15}, p_{16}, p_{17}, p_{18}, p_{19}, p_{20}.$$

Result is algebraically independent over $\mathbb{C}$: yes.

13. HIV [130, Page 1].

Original generating set information: 10 indeterminates; 77 non-constant functions; maximal total degrees of numerator and denominator are (8, 0); 2.3 KB total in string representation.

Indeterminates: $a, b, c, d, h, k, q, u, \lambda, \beta$.

Original generating set: too large to be listed.

Result of our algorithm:

$$a, b, d, h, u, cq^2, \frac{\lambda}{q}, \beta k q.$$

Result is algebraically independent over $\mathbb{C}$: yes.

14. HIV2 [114, Eq. (9)].

Original generating set information: 10 indeterminates; 9 non-constant functions; maximal total degrees of numerator and denominator are (4, 0); 221 bytes total in string representation.

Indeterminates: $b, c, d, s, k_1, k_2, q_1, q_2, w_1, w_2$.

Original generating set:

$$b, d, -s, -bk_2 q_2, b^2 k_2 q_2, c + k_1 + w_1 + w_2, -bk_2 q_2 s + ck_1 w_2 + cw_1 w_2, ck_1 + cw_1 + cw_2 + k_1 w_2 + w_1 w_2, bdk_2 q_2 - bk_1 k_2 q_1 - bk_1 k_2 q_2 - bk_2 q_2 w_1.$$

Result of our algorithm:

$$b, d, s, k_2 q_2, c + k_1 + w_1 + w_2, ck_1 w_2 + cw_1 w_2, k_1 k_2 q_1 + k_1 k_2 q_2 + k_2 q_2 w_1, ck_1 + cw_1 + cw_2 + k_1 w_2 + w_1 w_2.$$

Result is algebraically independent over $\mathbb{C}$: yes.
15. Influenza_MB1 [129, Figure 3 and Appendix A].

Original generating set information: 7 indeterminates; 27 non-constant functions; maximal total degrees of numerator and denominator are (4, 0); 403 bytes total in string representation.

Indeterminates: c, p, r, d<sub>I</sub>, k<sub>T</sub>, s<sub>T</sub>, β.

Original generating set:

$$
\begin{align*}
& c, d_I, k_T, k_T, -d_I, -k_T, -s_T, -\beta, ck_T, -ck_T, -\beta c, d_I k_T, -d_I k_T, -k_T s_T, -\beta d_I, -cd_I s_T, \\
& \frac{11}{1000} k_T, -d_I k_T s_T, -r + \frac{11}{1000}, -\beta k_T - c, -cd_I k_T s_T, \frac{11}{1000} d_I k_T, \frac{11}{1000} cd_I k_T, \\
& -\beta ck_T - d_I s_T, -\beta d_I k_T - d_I r + \frac{11}{1000} d_I, -\beta cd_I k_T - cd_I r + \frac{11}{1000} cd_I, -\beta cd_I.
\end{align*}
$$

Result of our algorithm:

$c, r, d_I, k_T, s_T, \beta.$

Result is algebraically independent over $\mathbb{C}$: yes.

16. Influenza_MB2 [129, Figure 3 and Appendix A].

Original generating set information: 8 indeterminates; 50 non-constant functions; maximal total degrees of numerator and denominator are (4, 0); 817 bytes total in string representation.

Indeterminates: c, p, r, c<sub>T</sub>, d<sub>I</sub>, k<sub>T</sub>, s<sub>T</sub>, β.

Original generating set: too large to be listed.

Result of our algorithm:

$c, r, c_T, d_I, k_T, s_T, \beta.$

Result is algebraically independent over $\mathbb{C}$: yes.

17. Influenza_MB3 [129, Figure 3 and Appendix A].

Original generating set information: 8 indeterminates; 27 non-constant functions; maximal total degrees of numerator and denominator are (4, 0); 372 bytes total in string representation.

Indeterminates: c, p, r, d<sub>I</sub>, d<sub>T</sub>, k<sub>T</sub>, s<sub>T</sub>, β.

Original generating set:

$$
\begin{align*}
& c, d_I, k_T, k_T, -d_I, -k_T, -s_T, -\beta, ck_T, -ck_T, -\beta c, d_I k_T, d_T - r, d_T k_T, -d_I k_T, \\
& -k_T s_T, -\beta d_I, -cd_I s_T, -\beta cd_I, d_I d_T k_T, -d_I k_T s_T, -\beta k_T - c, cd_I d_T k_T, -cd_I k_T s_T, \\
& -\beta ck_T - d_I s_T, -\beta d_I k_T + d_I d_T - d_I r, -\beta cd_I k_T + cd_I d_T - cd_I r.
\end{align*}
$$

Result of our algorithm:

$c, r, d_I, d_T, k_T, s_T, \beta.$

Result is algebraically independent over $\mathbb{C}$: yes.

18. Influenza_MB4 [129, Figure 3 and Appendix A].

Original generating set information: 9 indeterminates; 50 non-constant functions; maximal total degrees of numerator and denominator are (4, 0); 755 bytes total in string representation.

Indeterminates: c, p, r, c<sub>T</sub>, d<sub>I</sub>, d<sub>T</sub>, k<sub>T</sub>, s<sub>T</sub>, β.

Original generating set: too large to be listed.

Result of our algorithm:

$c, r, c_T, d_I, d_T, k_T, s_T, \beta.$

Result is algebraically independent over $\mathbb{C}$: yes.

19. Influenza_MD1 [129, Figure 3 and Appendix A].

Original generating set information: 6 indeterminates; 14 non-constant functions; maximal total degrees of numerator and denominator are (3, 0); 206 bytes total in string representation.

Indeterminates: c, p, r, d<sub>I</sub>, s<sub>T</sub>, β.

Original generating set:

$$
c, -c, d_I, -d_I, \beta, \beta c, d_I s_T, cd_I s_T, \frac{-1000}{11} r, \frac{-1000}{11} s_T, \frac{-11}{1000} d_I, \frac{-11}{1000} cd_I, \beta d_I + d_I r, \beta cd_I + cd_I r.
$$

Result of our algorithm:

$c, r, d_I, s_T, \beta.$

Result is algebraically independent over $\mathbb{C}$: yes.