\section{Introduction}
\smallskip
{\em Context and motivation.}
Enumerative combinatorics contains a vast landscape of non-trivial counting problems
that can hardly be solved without introducing the associated generating 
functions~\cite{Ardila15}. 
Setting up and solving a functional equation for such a generating function 
generally makes it possible to deduce properties of the discrete objects of interest 
(e.g. closed formulas~\cite[p.~$104$]{Drmota97},
asymptotic behavior~\cite[p.~147]{FlSe09}, \cite[\S 5.4.5]{Schaeffer15}).
Among those problems, many require refining the initial enumeration in order
to more easily write a functional equation~\cite{Brown64, Brown65a, Tutte68}. 
Algebraically, this process amounts to introducing
an additional variable, called~\emph{catalytic}~\cite{Zeilberger00}.\\
\indent The so-obtained functional equation
with one catalytic variable usually relates the refined generating function to
its specializations with respect to the catalytic variable, one of these specializations
being the generating function of the initial enumeration problem.
A standard way to deduce combinatorial properties from the functional equation is first to
determine the nature of the generating function~\cite[Chapter~$6$]{Stanley99}
(e.g. rational, algebraic, D-finite, \ldots),
then to compute a witness (e.g. an annihilating polynomial or an annihilating linear differential equation, \ldots).\\

\smallskip {\em A typical example.} One enumeration problem consists in studying bicolored maps
(black, white) such that the degree of each black face is~$3$ and the degree of each white
face is a multiple of~$3$. Such objects are called~\emph{$3$-constellations}~\cite[Fig.~$6$]{BMJ06};
they have been enumerated via a bijective approach in~\cite{BMS00}.
Let us define~$a_n$ as the number of~$3$-constellations with~$n$ black faces.
We consider the refined sequence~$a_{n, d}$ as the number of~$3$-constellations
having~$n$ black faces and outer degree~$3d$.
Let~$F(t, u) = \sum_{n, d\geq 0}{a_{n, d}u^dt^n}\in\mathbb{Q}[u][[t]]$ be its generating function\footnote{$F(t,u)$ has
polynomial coefficients in~$u$ since for a fixed number of black faces, the outer degree
is finite.}.
The catalytic variable is~$u$ and
the specialization~$F(t, 1)$ is the generating function of the sequence~$(a_n)_{n\geq 0}$.
Using a classical ``deletion of the root edge'' argument~\cite[Fig.~$7$]{BMJ06}, one gets the
following functional equation with 1 catalytic variable:
\begin{align}\label{eqn:funceqn_3const}
F(t, u) = 1 \;&+ tu(2F(t, u) + F(t, 1))\frac{F(t, u) - F(t, 1)}{u-1}\\
\nonumber  &+ tuF(t, u)^3 + \; tu\frac{F(t, u) - F(t, 1) - (u-1)\partial_uF(t, 1)}{(u-1)^2}.
\end{align}
Note that for any~$a\in\mathbb{Q}$, the \emph{divided difference
operator}~$\Delta_a:F\mapsto (F(t, u) - F(t, a))/(u-a)$ maps~$\mathbb{Q}[u][[t]]$ to itself.
As a consequence, it follows that Equation~\eqref{eqn:funceqn_3const} admits a unique
solution in~$\mathbb{Q}[u][[t]]$: the first fraction is $\Delta_1F$ while the second one is
$\Delta_1^2F\equiv \Delta_1(\Delta_1F)$.
Rewritten as~$F = 1 + tu(3F - (u-1)\Delta_1F)\Delta_1F + tuF^3 + tu\Delta_1^2F$,
Equation~\eqref{eqn:funceqn_3const} is a~\emph{discrete differential equation} (DDE)
of order~$2$, since the operator~$\Delta_1$ is iterated~$2$ times.\\
Equation~\eqref{eqn:funceqn_3const} has the
property that its unique solution~$F\in\mathbb{Q}[u][[t]]$ has a specialization~$F(t, 1)$
which is algebraic over~$\mathbb{Q}(t)$. More precisely, in~\eqref{eqn:funceqn_3const}
the specialization~$F(t, 1) = 1 + t + 6t^2 + 53t^3 + 594t^4 + \cdots$ is a root in~$x$
of~$-3t^2+47t-1+(27t^2-66t+1)x+(-81t^2+18t)x^2+81t^2x^3$.


\smallskip {\em A general algebraicity result.}
The algebraicity of the unique solutions in~$\mathbb{Q}[u][[t]]$
of~\eqref{eqn:funceqn_3const}
is in fact a consequence of the following strong and elegant result proved in~$2006$ by
Bousquet-Mélou and Jehanne. It ensures algebraicity of solutions
of the most frequent class of DDEs of arbitrary order~$k$ and with one catalytic variable, 
namely the class of DDEs of the fixed-point type.
  \begin{theorem} \emph{(\cite[Thm.~$3$]{BMJ06})}\label{thm:BMJ06}
Let~$\mathbb K$ be a field of characteristic~$0$
and consider two polynomials
$Q \in \mathbb K[x, y_1, \ldots, y_k, t, u]$
and
$f \in \mathbb K[u]$,
where $k\in\mathbb{N}\setminus \{ 0 \}$.
Let $a\in\mathbb K$ and
$\Delta_a :\mathbb K[u][[t]]\rightarrow \mathbb K[u][[t]]$ be the
divided difference operator $\Delta_a F(t, u) := (F(t, u) - F(t,
a))/(u-a)$.
Let us denote by $\Delta_a^{i}$ the operator obtained by applying $i$ times $\Delta_a$.
Then, there exists a unique solution $F\in \mathbb K[u][[t]]$ of the  equation
\begin{equation}\label{eqn:initfunceqn}
F(t, u) = f(u) + t \, Q(F(t, u), \Delta_a F(t, u), \ldots, \Delta_a^{k}F(t, u), t, u),
\end{equation}
and moreover $F(t,u)$ is algebraic over $\mathbb K(t, u)$.
  \end{theorem}
  %

This result~\cite[Thm.~$3$]{BMJ06} has been further extended in~\cite{NoSe22}
to the case 
of \emph{systems of DDEs} of the form~\eqref{eqn:initfunceqn} with 1
catalytic variable.\\
\indent In fact, the algebraicity results proved in~\cite{BMJ06, NoSe22} are particular cases
of a much deeper and older result in commutative algebra proved by Popescu~\cite{Popescu86} in
the context of Artin approximation theory with nested conditions.
The strength of the approaches presented in~\cite{BMJ06, NoSe22} 
lies in the effectivity of their algebraicity proofs. 
Despite a recent improvement in the linear case~\cite{CJPR19}, Popescu's result
is still not known to admit an effective proof.
 
  

\smallskip    \indent{\em Setting and main goal.}
    In the remainder of this article, we focus only on DDEs of the form~\eqref{eqn:initfunceqn}.
    Note that one can consider the associated \emph{polynomial functional equation} obtained
by multiplying \eqref{eqn:initfunceqn} by the smallest power of $(u-a)$ such that
the product is polynomial in~$u$. This new equation is denoted by
\begin{equation}\label{eqn:initpoleqn}
P(F(t, u), F(t, a), \ldots, \partial_u^{k-1}F(t, a), t, u) = 0,
\end{equation}
for some non-zero polynomial~$P\in\mathbb{K}[x, z_0, \ldots, z_{k-1}, t, u]$.
Starting from~\eqref{eqn:initpoleqn}, our main goal is to compute a
non-zero~$R\in\mathbb{K}[t, z_0]$ such that~$R(t, F(t, a)) = 0$. Remark that setting
$u=a$ in~\eqref{eqn:funceqn_3const} or differentiating~\eqref{eqn:funceqn_3const} 
with respect to $t$ or $u$ does not bring any non-trivial informations. 
In this article, we will focus on designing systematic algorithms for solving equations such as~\eqref{eqn:funceqn_3const}.

\smallskip {\em Previous work.}
We use the notation~$\overline{\mathbb{K}}[[t^{\frac{1}{\star}}]]\equiv
\cup_{d\geq 1}\overline{\mathbb{K}}[[t^{\frac{1}{d}}]]$.
For \emph{linear} DDEs, the \emph{kernel method}~ (which first appears
in an exercise of the famous book by Knuth~\cite[Ex.~2.2.1-4]{Knuth68}
and was systematized in~\cite{BaFl02})
consists in finding the roots
$u = U(t)\in\overline{\mathbb{K}}[[t^{\frac{1}{\star}}]]$ of the coefficient in~$x$ of~$P$.
By taking the resultant with respect to~$u$ of this coefficient and of~$P$, one obtains
a polynomial relation relating $F(t, a), \ldots, \partial_u^{k-1}F(t, a)$.
Since the work~\cite{BoPe00} of Bousquet-M\'elou and Petkov{\v{s}}ek, 
the linear case can be considered as fully understood.\\
An extension
of this method to the setting where~$\deg_x(P) = 2$ is also classical and is called
the~\emph{quadratic method}. It first appears in Brown's work~\cite{Brown65b}
for the case~$k=1$. It is based on a different elegant argument which
produces a polynomial relation between these specialized series.
This method was extended thirty years later on a particular
family of examples by Bender and Canfield~\cite{BeCa94}.\\
\indent
Both the kernel method and the quadratic method were generalized by the approach proposed
by Bousquet-Mélou and Jehanne in~\cite{BMJ06}.
Their method consists of starting from~\eqref{eqn:initpoleqn} and of creating more
polynomial equations having a non-trivial solution with~$F(t, a)$ as
its~$\{z_0\}$-coordinate. When this
method produces as many polynomials as variables and if the induced system generates
a $0$-dimensional ideal, a polynomial elimination strategy performed on the
polynomial system allows one to compute a non-zero element~$R\in\mathbb{K}[t, z_0]$
annihilating~$F(t, a)$. In the case where the system does not have the above properties,
a deformation of~\eqref{eqn:initfunceqn} via the introduction of a parameter
allows one to compute such an~$R$. \\
\indent This unified method contains however some intrinsic limitations due to the
number of variables introduced when creating more polynomial equations, and to the lack
of geometric interpretation of the problem. It was already mentioned in~\cite{BMJ06}
that the method was ``lacking an efficient elimination theory for polynomial systems
which (\ldots) are highly symmetric''.\\
A first step in the algorithmic study of DDEs has been initiated in~\cite{BoChNoSa22} for DDEs 
of order~$1$.
The case of DDEs of higher order was left as future work in~\cite{BoChNoSa22};
this is the purpose of the current article.

\smallskip {\em Main results.}
In~\cref{sec:generic_modeling}, we recall the polynomial system reduction from~\cite{BMJ06}
and provide a geometric interpretation of it.
In~\cref{sec:direct}, we prove under genericity assumptions on the input DDE~\eqref{eqn:initpoleqn}
that the algebraicity degree of~$F(t,a)$ is bounded
by~$\delta^{3k}/k!$~(Prop.~\ref{prop:bounds_under_H1}).
Here and in all that follows, $\delta$ denotes an upper bound on the total degree of~$P$  in~\eqref{eqn:initpoleqn}.
We deduce from this bound that a nonzero
annihilating polynomial of~$F(t, a)$ can be computed
in~$\tilde{O}(\delta^{6k}(k^2\delta^{k+3} + k^{k+1} + \delta^{1.89k}/k!))$ ops.
in~$\mathbb{K}$ (Prop. \ref{prop:complexity_under_H1}).
Here, and in the whole paper, the soft-O notation~$\tilde{O}(\cdot)$ 
hides polylogarithmic factors in the argument.
 In~\cref{sec:guessandprove},
we use the upper-bound $\delta^{3k}/k!$ to generalize Prop.~2.11 in~\cite{BoChNoSa22}
and deduce a complexity estimate
in~$\tilde{O}({k\cdot\delta^{10.12\cdot k}}/{(k-1)!^2})$ ops. in~$\mathbb{K}$
(Prop.~\ref{prop:hybridgp}).
In~\cref{sec:elimination}, we introduce a new
algorithm based on algebraic elimination and 
specialization properties of Gr\"obner bases. 
In~\cref{sec:stickelberger}, we
design one more new algorithm based on a geometric interpretation of the problem.
In~\cref{sec:experiments}, we provide
experimental results based on efficient implementations
of~\cref{sec:direct,sec:guessandprove,sec:elimination,sec:stickelberger}
for several DDEs coming from combinatorics.
The practical gains compared to the state-of-the-art go from a few minutes to several days,
 and we solve 
 the DDE related to~$5$-constellations ($k=4$)
 using a combination of~Sections~\ref{sec:guessandprove} and~\ref{sec:elimination}.

\smallskip{\em Notation.} 
We denote by~$\mathbb{K}$ an effective field of characteristic~$0$.
We write $\overline{\mathbb K}$ for an algebraic
closure of $\mathbb{K}$, and $\mathbb K[t], \mathbb K(t)$ and $\mathbb K[[t]]$ 
for, respectively, the rings of polynomials, rational functions and 
formal power series 
in~$t$ with coefficients in~$\mathbb K$.
We also use the notation
${\mathbb K}[[t^{1/\star}]]:=\bigcup_{d\geq1} \mathbb K[[t^{1/d}]]$ for the ring of~``fractional power
series'', that is series of the form 
$\sum_{n\geq 0} u_n t^{n/d}$ for some
integer $d\geq 1$.
We use the convention $\partial_x f$ and alike for the partial derivative of a function $f$
with respect to~$x$, we denote by $\operatorname{disc}_x(p)$
discriminant, with respect to a variable~$x$. 
For~$p$ a polynomial in~$n$ variables over~$\mathbb K$,
we denote 
by~$\deg(p)$ its total degree
and by~$\deg_x(p)$ its partial degree w.r.t the variable~$x$.
For an ideal $I \subset \mathbb K[x_1,\ldots, x_n]$, we denote by
$V(I)= V_{\overline{\mathbb K}}(I)$
the affine variety defined by $I$ over $\overline{\mathbb K}$.
 

\smallskip {\em Polynomial elimination basics.}  
We will repeatedly make use of the following 
fundamental results in polynomial elimination theory.
 For proofs and further context, we refer to 
 \cite[Chap.~3, Thm.~2, p.~122]{CoLiOS07} for Fact~\ref{fact:elim}\eqref{fact:elimination},
 \cite[Chap. 3, \S 5, Theorem 3, p. 159]{CoLiOS07} for Fact~\ref{fact:elim}\eqref{fact:extension},
\cite[Theorem~1.1]{Cox20} for Fact~\ref{fact:elim}\eqref{fact:eigenvalue}
and
\cite[Theorem~1.2]{Cox20} for Fact~\ref{fact:elim}\eqref{fact:Stickelberger}.

\begin{fact}\label{fact:elim}
Let $\mathbb K$ be a field, let $I \subset \mathbb K[x_1,\ldots, x_n]$ be an ideal.

\begin{enumerate}
	 \item (Elimination theorem) \label{fact:elimination}
Let $G$ be a Gr\"obner basis of $I$ with respect to the \emph{lex} order with $x_1 > \cdots > x_n$. Then, for every $0 \leq \ell \leq n$, the set
$G_\ell = G \cap \mathbb K[x_{\ell+1},\ldots,x_n]$ is a Gr\"obner basis of the $\ell$-th elimination ideal $I_\ell = I \cap \mathbb K[x_{\ell+1},\ldots,x_n]$.

\item (Extension theorem) \label{fact:extension}
Assume that $\mathbb K$ is algebraically closed.
Let $I=(f_1,\ldots,f_s)$, with  $f_i = c_i(x_2,\ldots,x_n)x_1^{N_i} + $ (lower degree terms in $x_1$). If $(a_2,\ldots,a_n) \in V(I_1)\setminus V(c_1,\ldots,c_s)$,
then there exists $a_1\in\mathbb K$ such that $(a_1,\ldots,a_n)\in V(I)$.
	
\item (Eigenvalue Theorem)\label{fact:eigenvalue}
Assume $I$ is zero-dimensional, let $\mathbb{A} = \mathbb K[x_1,\ldots, x_n]/I$
and $m_f : \mathbb{A} \rightarrow \mathbb{A}$ the multiplication-by-$f$ endomorphism of $\mathbb{A}$.
Then, 
the eigenvalues of $m_f$ are the values of $f$ at the finitely many points of $V(I)$.

\item (Stickelberger's theorem)\label{fact:Stickelberger}
If $I$ is radical and under the assumptions of \eqref{fact:eigenvalue},
the characteristic polynomial $\det(x I - m_f)$ of $m_f$ is equal to
$\prod_{a\in V(I)} (x-f(a))$.

\end{enumerate}
\end{fact}	

\smallskip {\em Complexity basics.}  
The algorithmic costs are estimated by counting elementary arithmetic operations
$(+, -, \times, \div)$ in the base field~$\mathbb{K}$ at unit cost. The
notation~$\theta$ refers to any feasible exponent for matrix multiplication over~$\mathbb K$.
The best current upper-bound is
$\theta<2.37188$~\cite{DuWuZh22}.
All classical operations on univariate polynomials of degree~$d$ in~$\mathbb{K}[x]$
(multiplication, multipoint evaluation and interpolation, extended gcd, resultant, squarefree part, etc) can be
performed in~$\tilde{O}(d)$.
We refer to the  book by von zur Gathen
and Gerhard~\cite{GaGe13} for these facts and related questions.



\section{From combinatorics to polynomials}\label{sec:generic_modeling}

\subsection{Solving DDEs via polynomial systems} 

The method of Bousquet-Mélou and Jehanne \cite[Section~$2$]{BMJ06} is based on the idea
of creating,
starting from the input equation~\eqref{eqn:initpoleqn}, 
new polynomial equations inducing a polynomial system
that admits a solution which has $F(t, a)$ as its $z_0$-coordinate.\\
The corresponding procedure is the following. 
One takes the derivative of \eqref{eqn:initpoleqn} with respect to the catalytic variable~$u$ and finds
\begin{align}\label{eqn:derivative_poleqn}
  \partial_uF(t, u)\;\cdot&\;
  \partial_xP(F(t, u), F(t, a), \ldots, \partial_u^{k-1}F(t, a), t, u)\\
\nonumber & + \partial_uP(F(t, u), F(t, a), \ldots, \partial_u^{k-1}F(t, a), t, u) = 0.
\end{align}

  Now,
  for any solution $u = U(t)$ in $\overline{\mathbb{K}}[[t^{\frac{1}{\star}}]]\setminus{\{a\}}$
  of the equation
\begin{equation}\label{eqn:crucial_eqn_BMJ}
\partial_xP(F(t, u), F(t, a), \ldots, \partial_u^{k-1}F(t, a), t, u) = 0,
\end{equation}
one obtains by plugging $u = U(t)$ in~\eqref{eqn:derivative_poleqn} that
\begin{align*}
  (x, z_0, \ldots, z_{k-1}, u) &=
  (F(t, U(t)), F(t, a), \ldots, \partial_u^{k-1}F(t, a), U(t))\\
  &\in\overline{\mathbb{K}}[[t^{\frac{1}{\star}}]]\times \mathbb{K}[[t]]\times \cdots
  \times \mathbb{K}[[t]]\times \overline{\mathbb{K}}[[t^{\frac{1}{\star}}]]
\end{align*}
is a solution of the following polynomial system
\begin{equation}\label{eqn:initial_system}
 \begin{cases} P(x, z_0, \ldots, z_{k-1}, t, u) = 0,\\
  \partial_xP(x, z_0, \ldots, z_{k-1}, t, u) = 0,\;\;\; u\neq a.\\
  \partial_uP(x, z_0, \ldots, z_{k-1}, t, u) = 0,\\
  \end{cases}
  \end{equation}

Showing the existence of solutions $u = U(t)$ in
$\overline{\mathbb{K}}[[t^{\frac{1}{\star}}]]\setminus{\{a\}}$ to~\eqref{eqn:crucial_eqn_BMJ}
  is generally done in combinatorics by applying
  \cite[Theorem~$2$]{BMJ06} and by computing the first terms of the solutions to show that
  they are not constantly equal to $a$.
  In order to avoid those checks, we make the following assumption.
\begin{align*}
\label{hyp1} \tag{\bf H1}
\underline{\textsf{Hypothesis 1:}} \quad
  &\deg_u(\partial_xP(x, z_0, \ldots, z_{k-1}, 0, u))\geq k,\\
  &\partial_{y_k}Q(f(a), f'(a), \ldots, \frac{1}{k!}f^{(k)}(a), 0, a)\neq 0.
\end{align*}


Under~\eqref{hyp1}, we show that there exist~$k$ distinct solutions~$u$
to the constraints given by~\eqref{eqn:initial_system}.
Note that~\eqref{hyp1} holds for a generic choice of $f$ and $Q$
in~\eqref{eqn:initfunceqn}.


\begin{prop}\label{prop:root_existence}

  Under assumption~\eqref{hyp1}, there exist $k$ distinct solutions
  $u = U_1(t), \ldots, U_k(t)$ of~\eqref{eqn:crucial_eqn_BMJ}.
  All of them lie in $\mathbb{L}[[t^\frac{1}{k}]]$,
  where $\mathbb{L}/\mathbb{K}$ is a field extension of degree upper bounded by $k$.

\end{prop}
\begin{proof}
  
    By expressing equation~\eqref{eqn:crucial_eqn_BMJ} in terms of the derivatives of~$Q$
    and by using the first part of~\eqref{hyp1},
    searching solutions $u \neq a$ of~\eqref{eqn:crucial_eqn_BMJ}
    is seen to be equivalent to looking for solutions $u\neq a$ of
 \begin{align*}
   (u&-a)^k
   =\,\, t\cdot(u-a)^k
   \cdot\partial_x Q(F(t, u), \Delta_a F(t, u), \ldots, \Delta_a^k F(t, u), t, u)\\
   \nonumber  &+
   t\cdot \sum\limits_{i = 1}^{k}{(u-a)^{k-i}\partial_{y_i}Q(F(t, u), \Delta_a F(t, u),
      \ldots, \Delta_a^k F(t, u), t, u)}.
 \end{align*}
 
   By specializing this equation
   at $t=0$, it follows that the
   constant term of any solution in $u$ of~\cref{eqn:crucial_eqn_BMJ} is equal to $a$.
   Then,
   taking the $k$-{th} root of the previous equation yields an identity of the form
   $u = a + t^{\frac{1}{k}}(\alpha + X)^{\frac{1}{k}}$, where
   $\alpha := \partial_{y_k}Q(f(a), f'(a), \ldots, \frac{1}{k!}f^{(k)}(a), 0, a)
   \in\mathbb{K}\setminus\{0\}$ and $X\in\mathbb{K}[[t]][[u]]$ satisfies $X(0, a) = 0$.
   Applying Taylor's multivariate formula in $t$ and $u-a$,
   one expands $(\alpha + X)^{\frac{1}{k}}$
   in~$\mathbb{L}[[t]][[u]]$, with~$\mathbb{L} = \mathbb{K}(\alpha^{\frac{1}{k}})$.
   By a fixed-point argument applied to $u = a + t^{\frac{1}{k}}(\alpha + X)^{\frac{1}{k}}$,
   the $k$-{th} roots of~$\alpha$ induce
   the existence of~$k$ distinct solutions
   $u = U_1(t), \ldots, U_k(t)\in\mathbb{L}[[t^{\frac{1}{k}}]]\setminus\{a\}$
   to equation~\eqref{eqn:crucial_eqn_BMJ}, all of them lying in~$\mathbb{L}[[t^\frac{1}{k}]]$.
 %
\end{proof}


The main idea of~\cite{BMJ06} is that the existence of the distinct solutions
$U_1(t), \ldots, U_k(t)$ induce distinct pairs
  $(F(t, U_i(t)), U_i(t))\in\overline{\mathbb{K}}[[t^{\frac{1}{k}}]]^2$
  for every $1\leq i \leq k$.
  Hence the point
  \begin{align}\label{eqn:point_I}
    (&x_1, u_1, \ldots, x_k, u_k, z_0, \ldots, z_{k-1})
    =\\
    \nonumber &(F(t, U_1(t)), U_1(t), \ldots, F(t, U_k(t)), U_k(t), F(t, a),
    \ldots, \partial_u^{k-1}F(t, a))
  \end{align}
  is a solution of the \emph{duplicated} polynomial system
  \begin{equation}\label{eqn:duplicated_system}
    \forall \;1\leq i \leq k\;, 
    \begin{cases} P(x_i, z_0, \ldots, z_{k-1}, t, u_i) = 0,\\
  \partial_xP(x_i, z_0, \ldots, z_{k-1}, t, u_i) = 0,\\
  \partial_uP(x_i, z_0, \ldots, z_{k-1}, t, u_i) = 0,\\
  \end{cases}
  \end{equation}
  defined by $3k$ equations in $3k$ variables over $\mathbb{K}(t)$.
  Now, to avoid the irrelevant solutions of~\eqref{eqn:duplicated_system},
  we restrict our attention to the solutions of~\eqref{eqn:duplicated_system}
  that are not solutions of
  $\prod_{i\neq j}({u_i - u_j})\cdot\prod_{i}({u_i - a}) = 0$; we define 
  $\textsf{diag}\in\mathbb{K}[u_1, \ldots, u_k]$ as the left-hand side of this equation.
  %
  \begin{notation}
  We write $\underline{x}$ (resp. $\underline{u}$ and $\underline{z}$) for the
  variables  $x_1, \ldots, x_k$ (resp. $u_1, \ldots, u_k$ and $z_0, \ldots, z_{k-1}$) 
  and  $\mathcal{I}$ for the ideal
  of~$\mathbb{K}(t)[\underline{x}, \underline{u}, \underline{z}]$ generated
  by~\eqref{eqn:duplicated_system}.
  \end{notation}
  %
  With the extra condition that
  $\textsf{diag}\neq 0$, those~$3k$ equations and variables generically define
  a $0$-dimensional ideal over~$\mathbb{K}(t)$ and hence induce 
  a finite set of solutions. For a later effective use of this finiteness, we introduce
  the following regularity assumption.
    %
\begin{align*}
\label{hyp2} \tag{\bf H2}
\underline{\textsf{Hypothesis 2:}} \quad
&\text{The ideal }\mathcal{I}_\infty := \mathcal{I}:\textsf{diag}^\infty\text{ in }
\mathbb{K}(t)[\underline{x}, \underline{u}, \underline{z}]\\
&\text{ is radical }\text{and has dimension}~0 \text{ over } \mathbb{K}(t).
\end{align*}

      \indent
      We now show that eliminating all variables in~$\mathcal{I}_\infty$ except~$t$
      and~$z_0$ yields a non-zero annihilating polynomial of~$F(t, a)$.
  
  \begin{prop}\label{prop:trivial_elimination_property}
    Under~\eqref{hyp1} and~\eqref{hyp2}, if
    $R\in\mathcal{I}_\infty\cap\mathbb{K}[t, z_0]\setminus\{0\}$
    then~$R(t, F(t, a)) = 0$.
  \end{prop}
  \begin{proof}
    Under~\eqref{hyp1}, we apply Proposition~\ref{prop:root_existence} to justify that
    the point given by~\eqref{eqn:point_I} lies in~$V(\mathcal{I})$. Now by definition of
    a saturated ideal, there exists $m\in\mathbb{N}$ such that
    $\textsf{diag}^m\cdot R\in\mathcal{I}$ writes as an algebraic expression
    in the polynomials involved in~\eqref{eqn:duplicated_system}.
    Specializing this expression to the point given by~\eqref{eqn:point_I} and using
    Proposition~\ref{prop:root_existence}, the point given by~\eqref{eqn:point_I}
    does not annihilate $\textsf{diag}$. Hence it annihilates~$R$.
    Finally,  the dimension property in~\eqref{hyp2} implies
    that $\mathcal{I}_\infty\cap\mathbb{K}[t, z_0]$
    is not reduced to~$\{0\}$.
  \end{proof}

  \subsection{Geometric interpretation} 
  We now introduce a more geometric interpretation of the fact that
  \eqref{eqn:point_I} is a solution of~\eqref{eqn:duplicated_system}.
   Recall that a subset of ${\acfield{t}}^{k}$
is said to be \emph{constructible} if
it is  a finite union of open subsets of a closed subset of~${\acfield{t}}^{k}$.
Typically,
  a zero set defined by polynomial equations and inequations is constructible.
  Hence denoting by~$\mathcal{X}\subset\acfield{t}^{k+2}$ the zero locus associated to the
  constraints~\eqref{eqn:initial_system}, we have that~$\mathcal{X}$ is a constructible set.
  We now define new geometric objects and assumptions for any constructible set
  $\mathscr{W}\subset {\acfield{t}}^{k+2}$, and deduce simple properties in the
  case~$\mathscr{W} = \mathcal{X}$.
  Let us define the
canonical projection $\pi : (x, u, \underline{z})\in\acfield{t}^{k+2}
\mapsto(\underline{z})\in\acfield{t}^k$ onto
the $\underline{z}$-coordinate space. In the whole paper, we assume that:
\begin{align*}
  \label{hyp:finite} \tag{\bf F}
  \text{The restriction of $\pi$ to $\mathscr{W}$ has finite fibers}. 
\end{align*}
For ${{\bm{\alpha}}}\in \acfield{t}^{k}$, we denote by $\mathscr{W}_{{{\bm{\alpha}}}}\subset
\acfield{t}^{2}\times \acfield{t}^{k}$ the fiber $\pi^{-1}({{\bm{\alpha}}})\cap \mathscr{W}$,
by $\card{u}{\mathscr{W}}{{{\bm{\alpha}}}}$ the number of $u$-coordinates of the points in
$\mathscr{W}_{{\bm{\alpha}}}$. We set $\cfiber{k}{u}{\mathscr{W}} := \{\bm{\alpha}\in{\acfield{t}}^{k} \mid \card{u}{\mathscr{W}}{{{\bm{\alpha}}}} \geq k\}$. 
      \begin{lemma}
		If $\mathscr{W}\subset {\acfield{t}}^{k+2}$ is constructible, then
		$\cfiber{k}{u}{\mathscr{W}}$
        is 
		also
        constructible. Moreover, under~\eqref{hyp1},
		$\cfiber{k}{u}{\mathcal{X}}$ is not empty.
        \end{lemma}
      \begin{proof}
        By duplicating $k$ times the
        variables~$x, u$ and eliminating all these duplicated variables,
        we obtain by Fact~\ref{fact:elim}\eqref{fact:elimination}
        that~$\cfiber{k}{u}{\mathscr{W}}$ is constructible. \\
         \indent  Under~\eqref{hyp1},~\cref{prop:root_existence} implies
        that the system~\eqref{eqn:initial_system} admits (at least)
        $k$ solutions in~$\overline{\mathbb{K}}[[t^\frac{1}{\star}]]$ with
        same~$(\underline{z})$-coordinates, and distinct $u$-coordinates. This proves
		that
         $\cfiber{k}{u}{\mathcal{X}}$ is not empty.
      \end{proof}

        The aim of the new
        algorithm that we will introduce in
        Sec.~\ref{sec:elimination}
	is to compute a finite set of polynomial equations and inequations
        in~$\mathbb{K}(t)[\underline{z}]$
        whose zero set in~${\acfield{t}}^k$
        is $\cfiber{k}{u}{\mathcal{X}}$.
		
        We now consider the projection map~$\pi_{\check{z_1}}:(x, u, \underline{z})
        \in\acfield{t}^{k+2}
        \mapsto (z_0, z_2, \ldots, z_{k-1})\in\acfield{t}^{k-1}$, and we assume in the rest
        of this paper that the following assumption holds:
        \begin{align*}
  \label{hyp:stick} \tag{$\mathbf{\check{F}}$}
  \text{The restriction of $\pi_{\check{z_1}}$ to $\mathscr{W}$ has finite fibers}. 
\end{align*}

        We now introduce a set~$\mathcal{S}_k(\check{z_1}, \mathscr{W})$
		that will provide a second geometric interpretation of our problem, 
		and will yield a second algorithm, given in~\cref{sec:stickelberger}. 
		We thus define
        \begin{align*}
          \mathcal{S}_k(\check{z_1}, \mathscr{W}) := \{\alpha =
          \;&(\alpha_0, \ldots, \alpha_{k-1})\in\overline{\mathbb{K}(t)}^k|
          \;\; \alpha\in\pi(\mathscr{W}) \;\;\wedge\;\; \\
          &\#\mathscr{W} \cap  \pi_{\check{z_1}}^{-1}((\alpha_0, \alpha_2, \ldots, \alpha_{k-1})) \geq k\}.
        \end{align*}
        \begin{lemma}\label{lemma:SW_constructible}
          The set $\mathcal{S}_k(\check{z_1}, \mathscr{W})$ is 
          constructible. Moreover, under assumption~\eqref{hyp1}
		  the set $\mathcal{S}_k(\check{z_1}, \mathcal{X})$
          is not empty.
        \end{lemma}
        \begin{proof}
          Considering polynomial constraints defining the set~$\mathscr{W}$,
          the cardinality condition can be modeled by: fixing the variables $z_0, z_2, \ldots,
          z_{k-1}$, duplicating the variables~$x, u, z_1$ and defining appropriate
          inequations in order to ensure the solutions of such a system to be distinct.
           If~\eqref{hyp1} holds, then~$\mathcal{S}_k(\check{z_1}, \mathcal{X})$ is not empty since
          it contains the point~$(F(t, a), \ldots, \partial_u^{k-1}F(t, a))\in\acfield{t}^k$.
        \end{proof}
        
        In Sec.~\ref{sec:stickelberger}, we will introduce a new algorithm that computes
        a finite set of polynomials in~$\mathbb{K}(t)[\underline{z}]$ whose zero
        set is~$\overline{\mathcal{S}_k(\check{z_1}, \mathcal{X})}$.
        Note that~$\cfiber{k}{u}{\mathcal{X}}$
        and~$\mathcal{S}_k(\check{z_1}, \mathcal{X})$ are related as follows:

        \begin{lemma}
          The following inclusion holds~$\cfiber{k}{u}{\mathcal{X}}\subset
          \mathcal{S}_k(\check{z_1}, \mathcal{X})$.
        \end{lemma}
        \begin{proof}
          Let us choose~$\alpha = (\alpha_0, \ldots, \alpha_{k-1})\in
          \cfiber{k}{u}{\mathcal{X}}$.
          By definition of~$\cfiber{k}{u}{\mathcal{X}}$ we have
          that~$\alpha\in\pi(\mathcal{X})$, hence
          $(\alpha_0, \alpha_2, \ldots, \alpha_{k-1})\in\pi_{\check{z_1}}(\mathcal{X})$. Now,
          any of the~$k$ points in~$\mathcal{X}\cap\pi^{-1}(\alpha)$ also belongs to
          $\mathcal{X}\in\pi_{\check{z_1}}((\alpha_0, \alpha_2, \ldots, \alpha_{k-1}))$.
          Hence~$\alpha\in\mathcal{S}_k(\check{z_1}, \mathcal{X})$.
          \end{proof}

 \section{Direct approach: degree bounds and complexity}\label{sec:direct}

  %
 In this section, we focus on the complexity of computing a precise
 nonzero element $R$ of $\mathcal{I}_\infty\cap\mathbb{K}[t, z_0]$, using the
 work of Bousquet-Mélou and Jehanne~\cite{BMJ06}. The
 following analysis is a generalization of the ones
 given in~\cite[Proposition~$2.8$]{BoChNoSa22}. It takes benefit of the group
 action of the symmetric group~$\mathfrak{S}_k$ of degree $k$ on the
 zero set~$V(\mathcal{I}_\infty)$, which can be exploited for DDEs of
 order~$k>1$.
  %
  \begin{prop}\label{prop:bounds_under_H1}
  Let $P$ be as in \eqref{eqn:initpoleqn} of total degree $\delta$.
  Assume that~\eqref{hyp1} and \eqref{hyp2} hold.
  Then $\deg(\mathcal{I}_\infty)$ and $\#V(\mathcal{I_\infty})$ are bounded
  by~$\delta^k\cdot(\delta - 1)^{2k}$, and there exists
  $R\in\mathbb K[t, z_0] \setminus \{ 0 \}$
  satisfying $R(t, F(t, a))~=~0$, with $\deg_t(R)$ and $\deg_{z_0}(R)$ at most
  $\mathfrak{b}$, where
   $ \mathfrak{b}~:=~{\delta^k(\delta-1)^{2k}}/{k!}$.
\end{prop}
\begin{proof}
  First, we
  identify a nonzero polynomial in ${\mathcal{I}}_\infty\cap\mathbb K[t, z_0]$.
  Note that~\eqref{hyp1} and~\eqref{hyp2} allow us to
  apply~\cref{prop:trivial_elimination_property} which implies that such a polynomial
  annihilates~$F(t, a)$.
  Now,~\eqref{hyp2} implies that the quotient ring $\mathbb K(t)[\underline{x},
  \underline{u}, \underline{z}]/{\mathcal{I}}_\infty$ is a finite dimensional
  $\mathbb K(t)$-vector space. Hence by
  Fact~\ref{fact:elim}\eqref{fact:eigenvalue},
  the endomorphism $m_{z_0}:f\mapsto z_0\cdot f$ admits a
  characteristic polynomial $\xi_{z_0}\in\mathbb K(t)[z_0]$ whose roots are exactly
  the $z_0$-coordinates of all points in $V({\mathcal{I}}_\infty)$ (in finite
  number by assumption~\eqref{hyp2}). Hence, multiplying $\xi_{z_0}$ by the lcm of
  the denominators of
  its coefficients and denoting by~$R$ the squarefree part of the resulting polynomial,
  the radicality of $\mathcal{I}_\infty$ implies that~$R\in\mathcal{I}_\infty\cap
  \mathbb{K}[t, z_0]$. Hence~$R$ satisfies~$R(t, F(t, a)) = 0$.\\
  %
  \indent We prove now that the partial degrees of $R$ in $t$ and $z_0$ are both
  bounded by $\mathfrak{b}$. We apply the exact same proof as the one done
  for proving~\cite[Proposition~$2.8$]{BoChNoSa22} but with
  $\mathcal{I}$ replaced by $\mathcal{I}_\infty$ and with
  $\deg(\partial_xP)$ and~$\deg(\partial_uP)$ both bounded by~$\delta-1$. This implies
  that~$\deg(\mathcal{I}_\infty)$ and $\# V(\mathcal{I}_\infty)$ are bounded
  by $\delta^k(\delta-1)^{2k}$.
  As the partial degrees of~$R$ are bounded
  by~$\deg(\mathcal{I}_\infty)$, it follows that~$\deg_{z_0}(R)\leq \delta^k(\delta-1)^{2k}$.\\
  %
  \indent It remains to justify the division by $k!$.
  We use the group action of~$\mathfrak{S}_k$ on
  $V(\mathcal{I}_\infty)$ obtained by permuting the pairs of
  coordinates given by $(x_i, u_i)$ and $(x_j, u_j)$, for $i\neq j$.
  This group action stabilizes the $(z_0, \ldots, z_{k-1})$ coordinate subspace.
  Hence one spares the cardinality of the orbits of this group action
  by dividing $\deg(\mathcal{I}_\infty)$
  by~$k!$.
  Thus $\deg_{z_0}(R)$ is bounded by~$\mathfrak{b}$. By permuting the roles of~$z_0$ and~$t$,
  one also obtains that~$\deg_t(R)\leq \mathfrak{b}$.
  %
  \end{proof}

 
\begin{prop}\label{prop:complexity_under_H1}
  Let $P$ be as in \eqref{eqn:initpoleqn} of total degree $\delta$. We suppose
  that~~\eqref{hyp1} and~~\eqref{hyp2} hold. Then there exists an algorithm
  which takes as input a straight-line program of length $L$
  evaluating~$P$ and~$\textsf{diag}$ and returns a nonzero polynomial $R\in\mathbb K[t,
    z_0]$ 
  such that~$R(t, F(t,
  a))=0$,  using
  \begin{align*}
    \tilde{O}((kL+1)\delta^{2k}(\delta-1)^{4k} &+ \delta^{2.63k}(\delta-1)^{5.26k}/k!)\\
    &\subset
    \tilde{O}(\delta^{6k}(k^2\delta^{k+3} + k^{k+1} + \delta^{1.89k}/k!))
  \end{align*}
  arithmetic operations in~$\mathbb K$.
\end{prop}
\begin{proof}
  %
  We generalize the proof of~\cite[Proposition~$2.9$]{BoChNoSa22} to our
  situation. By~\cref{prop:bounds_under_H1},
  $\deg(\mathcal{I}_\infty)$ and $\#V(\mathcal{I_\infty})$ are bounded
  by~$\delta^k\cdot(\delta - 1)^{2k}$. Using the algorithm
  underlying~\cite[Theorem~$2$]{Schost03}, we thus compute
  a parametric geometric resolution \cite{Schost03} of the zero
  set~$V(\mathcal{I}_\infty)$
  in~$\tilde{O}((kL+1)\delta^{2k}(\delta-1)^{4k})$ ops. in~$\mathbb{K}$.
  This algorithm computes two polynomials~$V(t, \lambda), W(t, \lambda)\in
  \mathbb{K}(t)[\lambda]$ giving the
  parametrization~$z_0=V(t, \lambda)/\partial_\lambda W(t, \lambda)$ whenever
  $W(t, \lambda)=0$.
  Performing evaluation-interpolation on~$t$ with~$O(\mathfrak{b})$ points, 
  we then
  compute the squarefree part~$R$ of the
  minimal polynomial~$\xi_{z_0}$ of the map~$m_{z_0}:f\mapsto f\cdot z_0$
  in~$\mathbb{K}(t)[\underline{x}, \underline{u}, \underline{z}]/\mathcal{I}_\infty$
  by taking the bivariate resultant of~$z_0\cdot\partial_\lambda W(t, \lambda) -
  V(t, \lambda)$ and
  $W(t, \lambda)=0$. Using~\cite[\S 5]{HyNeSc19},
  this is done 
  in~$\tilde{O}(\delta^{2.63k}(\delta-1)^{5.26k}/k!)$
  ops.~in~$\mathbb{K}$. 

  Finally, the inclusion comes from the cost for evaluating~$P$
  and the saturating polynomial~$\textsf{diag}$, which by the Baur-Strassen
  theorem~\cite[Theorem~$1$]{BaSt83} satisfies~$L\in O(k\delta^{k+3} + k^k)$.
\end{proof}

Despite the process of duplicating variables as done for obtaining the
system~\eqref{eqn:duplicated_system} is, up to deforming the initial DDE,
very fruitful for creating zero dimensional ideals and showing theoretical
algebraicity results~\cite{BMJ06,NoSe22}, it usually suffers from efficiency issues.
\noindent
Applying~\cite[Prop.~$2$]{Heintz83} for~$n, m\in\mathbb{N}$ and an algebraic
set~$V\subset\overline{\mathbb{K}}^n$,
we have that~$\deg(V^m) = \deg(V)^m$. The degenerate behavior of the state of the art
when $k$ grows up comes from this exponential growth of the ideal's degree in the number of
duplications (which is $k$ in our case).
Moreover, duplicating variables also introduces $3k+2$ variables, while the initial system
in~\eqref{eqn:initial_system} only deals with $k+3$ variables.
A natural hope is hence to
avoid these duplications by
a careful analysis of the geometry given by the initial
constraints~\eqref{eqn:initial_system}.

\begin{remark}
  This number of duplications and the group action of~$\mathfrak{S}_k$
over~$V(\mathcal{I}_\infty)$
is usually exploited by the state-of-the-art of the
polynomial system solving theory by working
in the invariant ring associated to this group action (see~\cite{FaSv12}). However in our
case, this approach would imply to introduce a number of variables which would be
at least equal to~$k$.
In the algorithm we introduce in~\cref{sec:elimination,sec:stickelberger}, we focus
on introducing no more extra variables.
\end{remark}


\section{Hybrid guess-and-prove algorithm}\label{sec:guessandprove}




  \indent
  We analyze in what follows the complexity of the \emph{hybrid guess-and-prove} algorithm
introduced in \cite[\S2.2.2]{BoChNoSa22}.
Recall that it blends algebraic elimination with a guess-and-prove approach
inspired by Zeilberger's method \cite{Zeilberger92}, see \cite[\S2.2.1]{BoChNoSa22}.
For functional equations of arbitrary order, the motivation of this
method comes from certain concrete examples for which the involved polynomial systems are
difficult to solve (e.g. \cite[\S3.6]{BMWa20}). 
Let us recall the algorithm:


\smallskip \underline{Hybrid guess-and-prove method}
\begin{enumerate}
  \item[(0)] Compute $F(t,a)$ mod~$t^\sigma$ for some integer~$\sigma\geq 1$;
  \item[(1)] Guess $R \in \mathbb K[t, z_0]\!\setminus\!\{ 0 \}$ such that
    $R(t, F(t,a))=0 \bmod t^\sigma$;
  \item[(2)] Check if $R(t, F(t,a)) = 0 \bmod t^{\mathfrak{b} \cdot \deg R + 1}$;
  \item[(3)] If not, then go back to~(0) with $\sigma:=2\sigma$; if yes, then return~$R$.
\end{enumerate}
  The correctness of this method is a consequence of~\cite[\S2.2.2]{BoChNoSa22}
  and of the existence, under suitable hypothesis, of a nonzero
  polynomial $R\in\mathbb{K}[t, z_0]$
  annihilating $F(t, a)$ with partial degrees bounded by~$\mathfrak{b}$.
  It remains to make those hypotheses completely explicit.
  \begin{notation}
    We still assume the existence
  of~$k$ distinct non-constant solutions $u=U_1(t), \ldots, U_k(t)$
  to~\eqref{eqn:crucial_eqn_BMJ} and
  we denote by $\mathcal{P}$ the point of $\overline{\mathbb{K}}[[t^{\frac{1}{\star}}]]^{3k}$
  obtained by concatenating the values
  $\{F(t, U_i(t))\}_{1\leq i\leq k}, \{U_i(t)\}_{1\leq i\leq k}$
  and $\partial_u^iF(t, a)\}_{0\leq i\leq k-1}$. 
  \end{notation}
\vspace{-0.5cm}
\begin{align*}
\label{hyp4} \tag{\bf H4}
\underline{\textsf{Hypothesis 4:}} \quad
&\text{The Jacobian matrix } \textsf{Jac} \text{ of }
\eqref{eqn:duplicated_system},\\
&\text{ considered in } \underline{x},
\underline{u}, \underline{z}, 
\text{ is invertible at } \mathcal{P}.
\end{align*}
\begin{lemma}\label{lemma:bounds_under_H3}
  Under~\eqref{hyp1} and \eqref{hyp4}, the saturation
  $\mathcal{I}:\textsf{Jac}^\infty$ is a radical and $0$-dimensional ideal
  of $\mathbb{K}(t)[\underline{x}, \underline{u}, \underline{z}]$. Hence there exists
  a nonzero
  $R\in\mathcal{I}:\textsf{Jac}^\infty\cap\mathbb{K}[t, z_0]$
  annihilating $F(t, a)$ with partial degrees bounded by~$\mathfrak{b}$.
\end{lemma}
\begin{proof}
  The radicality and dimension results are consequences of~\cite[Lemma~$2.10$]{BoChNoSa22}.
  Now, the rest of the proof is the same as the one used for proving
  Proposition~\ref{prop:bounds_under_H1}, with $\mathcal{I}_\infty$ replaced by
  $\mathcal{I}:\textsf{Jac}^\infty$.
  \end{proof}
This concludes the correctness of the hybrid guess-and-prove method in the case of DDEs of
order~$k>1$. Under~\eqref{hyp1} and~\eqref{hyp4}, we also deduce a complexity estimate
generalizing~\cite[Prop.~2.11]{BoChNoSa22}.


\begin{prop}\label{prop:hybridgp}
 Let $P$ be as in~\eqref{eqn:initpoleqn} and let $\delta$~be its total degree.
 Let us assume that assumptions \eqref{hyp1} and~\eqref{hyp4} hold, and that there exists a
 straight-line program of length~$L$ evaluating~$P$.
Then the hybrid guess-and-prove method terminates on input~$P$ using
  \begin{equation*}
    \tilde{O}(L\cdot k^3\cdot \mathfrak{b}^2 +\,k\cdot \mathfrak{b}^{\theta+1})
    \subset
    \tilde{O}({k\cdot\delta^{10.12\cdot k}}/{(k-1)!^2})
  \end{equation*}
  arithmetic operations in $\mathbb K$.
\end{prop}
\begin{proof}
  %
  We analyze the last execution of steps $(0)$--$(3)$, happening with
  $\sigma = O(\mathfrak{b}^2)$ and $\deg(R) = O(\mathfrak{b})$.
  %
  Using Hermite-Padé approximation (e.g.~\cite{GiJeVi03,BoJeSc08}), the guessing at
  step $(1)$ is done in $\tilde{O}(\sigma\cdot\deg_{z_0}(R)^{\theta-1})
  \subset\tilde{O}(\mathfrak{b}^{\theta+1})$ ops. in~$\mathbb K$.
  The order of $F(t, a)$ and the truncation order $\mathfrak{b}\cdot\deg(R) + 1$ at
  step $(2)$ are in~$O(\mathfrak{b}^2)$. Hence the truncated evaluation of $R$ at $z_0 = F(t, a)$
  is computed in $\tilde{O}(\deg(R)\cdot(\mathfrak{b}\cdot\deg(R)))\subset
  \tilde{O}(\mathfrak{b}^3)\subset \tilde{O}(\mathfrak{b}^{\theta+1})$ ops. in~$\mathbb K$.
  
  
  We compute the truncated series $F(t, a)$ using the classical Newton method, by considering
  the iteration
  \begin{equation}\label{eq:newton}
\mathbf F \mapsto \mathcal{N}(\mathbf F) := \mathbf F - J(\mathbf F)^{-1} \cdot S^{\mathrm{T}} (\mathbf F) ,
\end{equation}
  where $V$ is the system given by \eqref{eqn:duplicated_system},
  $J = \textsf{Jac}(V)$ (w.r.t. $\underline{x}, \underline{u}, \underline{z}$)
and $\mathbf F \in \mathbb{L}[[t^\frac{1}{k}]]^{3k}$ denotes an approximation of the solution
$\mathcal{F}$ at which $J$ is invertible by assumptions~\eqref{hyp1} and~\eqref{hyp4},
with $\mathbb{L}/\mathbb{K}$ a finite field extension
of degree at most~$k$ (by Proposition~\ref{prop:root_existence}).


Before going further, recall that the cost for elementary arithmetic operations in $\mathbb{L}$
can be expressed in terms of elementary arithmetic operations in~$\mathbb{K}$: as the field extension is
of degree at most~$k$, multiplying and summing two elements of~$\mathbb{L}$ can be done in
$\tilde{O}(k)$ arithmetic operations in~$\mathbb{K}$.


As the number of correct terms is doubled at each iteration
of~\eqref{eq:newton}, we iterate \eqref{eq:newton} at most
$O(\log(k))$ times to obtain a truncation of 
$\mathcal{F} \bmod t$ (because of the ramification appearing in the Puiseux series~$U_i$).
Hence we perform this precomputation before doubling the integer
power in $t$.


By the Baur-Strassen theorem \cite{BaSt83}, a straight-line program of length $O(L\cdot k)$
evaluating $V$ can be obtained from the one evaluating~$P$. By iterating this argument,
one also finds a straight-line program of length $O(L\cdot k)$
that evaluates the Jacobian matrix.
Consequently, there is also one in $O(L\cdot k)$ for its inverse.
Evaluating~$V$ and $J^{-1}$ at~$\mathbf F$ of some order $N$ requires consequently
$\tilde{O}(L\cdot k\cdot N)$ arithmetic operations in $\mathbb{L}$,
because of the fact that the
cardinality of the support of
$U_i \bmod t^N$ is in $O(k\cdot N)$. This is also the cost of a Newton iteration.
All in all, one obtains a complexity for step $(0)$ in $\tilde{O}(L\cdot k^2\cdot \sigma)$
arithmetic operations in $\mathbb{L}$.
Summing all complexities, the hybrid guess-and-prove complexity is
in~$\tilde{O}(L\cdot k^2\cdot \mathfrak{b}^2 +\, \mathfrak{b}^{\theta+1})$ ops. in $\mathbb{L}$.
This gives the global complexity
of~$\tilde{O}(L\cdot k^3\cdot \mathfrak{b}^2 +\, k\cdot \mathfrak{b}^{\theta+1})$ ops. in~$\mathbb{K}$.
The inclusion is a consequence of the estimate $L\in O(\delta^{k+3})$.
\end{proof}



\section{Approach using elimination theory}\label{sec:elimination}


\def\mathcal{H}{\mathcal{H}}
\def\mathscr{W}{\mathscr{W}}

Let $\mathscr{W}\subset \overline{\mathbb{K}(t)}^{k+2}$ be a 
{constructible} set defined by polynomials in $\mathbb{K}(t)[x, u, z_0, \ldots, z_{k-1}]$.
We assume that the finite fibers assumption~\eqref{hyp:finite} holds.
For $i\in \mathbb{N}$, we consider the (possibly infinite) set
\[
\cfiber{i}{u}{\mathscr{W}} \coloneqq
  \{\bm{\alpha}\in{\acfield{t}}^{k} \mid \card{u}{\mathscr{W}}{{{\bm{\alpha}}}} \geq i\}.
\]

Observe that $\cfiber{i + 1}{u}{\mathscr{W}}\subseteq \cfiber{i}{u}{\mathscr{W}}$ for any
$i\geq 1$.

Recall that the zero locus of a set of equations and inequations is a constructible set.
Thus one deduces that $\cfiber{i}{u}{\mathscr{W}}$ is a constructible set
by duplicating the variables
$x, u$ and the polynomial constraints defining~$\mathscr{W}$ in order to ensure a given cardinality
for the fiber of~$\pi$.
This is proved by combining the Rabinowitsch trick
with Fact~\ref{fact:elim}\eqref{fact:elimination}.

In this section, we provide an algorithm that takes as input $i\in
\mathbb{N}$ and a polynomial system defining some algebraic set $W\subset
\acfield{t}^{k+2}$ and that computes a set of polynomial equations and
inequations in $\field{t}[\underline{z}]$ whose solution set in
$\acfield{t}^{k}$ is~$\cfiber{i}{u}{W}$.
We then show how to apply this algorithm to compute witnesses of
algebraicity to solutions of DDEs of order~$k$. 


To begin with, we assume that $W$ is given by a polynomial system ${{\bm{f}}}$ in
$\mathbb{K}[t][x, u, \underline{z}]$ and we denote by $\mathcal{J}$ the ideal it
generates in~$\mathbb{K}(t)[x, u, \underline{z}]$.
To design our algorithm, we leverage advanced results of the theory of Gröbner
bases to characterize $\cfiber{i}{u}{W}$.

We start with a transfer lemma. We refer to \cite{CoLiOS07} for preliminaries on
monomial orderings and Gr\"obner bases.


Let $\rho_x: (x, u, \underline{z}) \mapsto (u,
\underline{z})$ be the canonical projection which forgets the variable~$x$. We
denote by $G$ a Gröbner basis for {$(\mathcal{J}, \succ)$}, 
where $\succ$ is a
lexicographic monomial ordering with $x \succ u \succ \uline{z}$. Let $G_x =
G\cap \field{t}[u, \uline{z}]$ and $\ell_x$ be the leading coefficients w.r.t.
the variable $x$ of the polynomials in $G$ which have positive partial degree
w.r.t. $x$. Finally, let $\mathscr{W}$ be $\rho_x(W)$.

\begin{lemma}\label{lemma:transfer}
  The set $\cfiber{i}{u}{W}$ coincides with $\cfiber{i}{u}{\mathscr{W}}$. 
\end{lemma}

\begin{proof}
  Let ${{\bm{\alpha}}}\in\pi(W)$. Observe that the $u$-coordinates of $W_{{\bm{\alpha}}}$
  coincide with the $u$-coordinates of the points in $\mathscr{W}$ which project on~${{\bm{\alpha}}}$. 
  Hence, the conclusion follows from the definitions of
  $\cfiber{i}{u}{W}$ and $\cfiber{i}{u}{\mathscr{W}}$.
\end{proof}

\begin{lemma}\label{lemma:proj_x}
  $\mathscr{W} = \rho_x(W)$ is defined by the vanishing of all polynomials in $G_x$ and
  the non-vanishing of at least one element in~$\ell_x$.
\end{lemma}

\begin{proof}
  This lemma comes by applying directly Fact~\ref{fact:elim}\eqref{fact:elimination} and
  Fact~\ref{fact:elim}\eqref{fact:extension}.
\end{proof}

We will use Lemmas~\ref{lemma:transfer} and~\ref{lemma:proj_x} to compute a
polynomial system that encodes $\cfiber{i}{u}{\mathscr{W}}$. Remark that by
applying Lemma~\ref{lemma:proj_x}, $\mathscr{W}$ is the union of the locally closed
sets defined by the vanishing of all polynomials in $G_x$ and the non-vanishing
of at least one element of $\ell_x$. Note that the vanishing set of $G_x$ is the
Zariski closure of $\mathscr{W}$. Furthermore, we denote by $G_u$ the set $G_x\cap
\field{t}[\uline{z}]$.

Given a set $F$ of polynomials in $\field{t}[u, \underline{z}]$ and $r\in
\mathbb{N}$, we denote by ${\sf DEG}_r(F, u)$ (resp. ${\sf DEG}_{\leq r}(F, u)$)
the subset of polynomials in $F$ of degree $r$ (resp. at most $r$) in $u$. For a
polynomial $f \in \field{t}[u, \underline{z}]$, we denote by ${\sf coeffs}(f,
u)$ its coefficients when $f$ is seen in $\mathbb{K}'[u]$ with
$\mathbb{K}' := \field{t}[\underline{z}]$.

\begin{lemma}
  We reuse the notation introduced above. If there is no polynomial in $G_x$
  whose partial degree w.r.t. $u$ is greater than or equal to $i$, then
  $\cfiber{i}{u}{\mathscr{W}}$ is empty.
\end{lemma}

\begin{proof}
  Suppose that~$\mathcal{F}_i(u, \mathscr{W}) \neq \emptyset$ and
  pick~$\alpha\in\mathcal{F}_i(u, \mathscr{W})$. 
  Denote
  by~$g\in G_u$ a non-zero polynomial of partial degree~$\leq i-1$ in~$u$.
  Suppose that there exists $1\leq j \leq \deg_u(g)$ such that the coefficient of~$u^j$
  does not vanish at~$\alpha$, and such that the coefficients of~$u^\ell$ vanish
  at~$\alpha$, for all~$\ell>j$. Then, there exists at most $j<i$ solutions to the
  equation~$g(u, \underline{z}=\alpha)=0$. This contradicts the fact that the fiber
  above~$\alpha$ has cardinality at least~$i$, and hence 
  that~$\alpha\in\mathcal{F}_i(u, \mathscr{W})$. 
\end{proof}

Let $\ell^{(i)}_u$ be the set of leading coefficients of the
polynomials in~$G_x$ that have partial degree at least $i$ w.r.t. $u$; we
denote $\ell_u^{(1)}$ by $\ell_u$.

For $g\in G_x$ of partial degree $i+j$ w.r.t. $u$, we denote by $c(g)$ the
Gröbner basis of the ideal $\langle g, \frac{\partial g}{\partial u}, \ldots,
\frac{\partial^{({j})} g}{\partial u^{j}} \rangle\cap \field{t}[\uline{z}]$.
Note that $g(u, {{\bm{\alpha}}})$ has at least $i$ distinct roots when ${{\bm{\alpha}}}$ does
not lie in the zero locus of $c(g)$.


Let $\mathcal{S}^{(i)}$ be the set of points $\mathbf{\beta}~\in\acfield{t}^{k+1}$ such that
the following polynomial constraints are simultaneously satisfied:
\begin{itemize}
\item[\em (a)] all polynomials in $G_u$ vanish at $\mathbf{\beta}$;
\item[\em (b)] all polynomials in ${\sf coeffs}(f, u)$ for $f\in {\sf DEG}_{\leq
    i-1}(G_u, u)$ vanish at $\mathbf{\beta}$;
\item[\em (c)] at least one resultant (w.r.t. $x$) of $g$ and $\ell$ does not vanish at~$\mathbf{\beta}$,
where $\ell\in \ell_x$ and $g$ is a
  polynomial of $G$ of positive partial degree w.r.t. $x$;
\item[\em (d)] at least one polynomial $\ell\in \ell^{(i)}_u$ does not vanish at~$\mathbf{\beta}$;
\item[\em (e)] at least one polynomial in $c(g)$ does not vanish at~$\mathbf{\beta}$,
for some $g \in G_x$ with
  partial degree at least $i$ w.r.t. $u$.
\end{itemize}

Further, we denote by {\em (f)} the disjunction $\vee_{\ell\in \ell_u} \ell\neq
0$. 

\begin{proposition}\label{prop:elim}
  The projection of~$\mathcal{S}^{(i)}$ onto the $\underline{z}$-coordinate space
  is~$\cfiber{i}{u}{\mathscr{W}}$.
\end{proposition}

\begin{proof}
  We start by showing that the set of points which do satisfy {\em (a)}, {\em
    (c)} and {\em (f)} coincides with $\pi(W)$ which, by definition, contains
  $\cfiber{i}{u}{W}$.

  Fact~\ref{fact:elim}\eqref{fact:elimination} shows that
  the Zariski closure of $\pi(W)$ is defined by~{\em~(a)}. 

  By Lemma~\ref{lemma:proj_x}, $\mathscr{W}$ is defined by the vanishing of all
  polynomials in $G_x$ and the non-vanishing of at least one element in
  $\ell_x$. Note that $\pi(W)$ coincides with the projection of $\mathscr{W}$ on the
  $\uline{z}$-coordinate space (which we also denote by $\pi(\mathscr{W})$ by a slight
  abuse of notation). Note that all points of $\mathscr{W}$ also satisfy {\em (c)}
  since those resultants would just vanish at the set of common solutions to
  $\ell, g$ where $\ell\in \ell_x$ and $g$ is a polynomial of $G$ of positive
  partial degree w.r.t. $x$. Hence their projection on the $\uline{z}$-space
  does also satisfy {\em (c)}. Applying
  Fact~\ref{fact:elim}\eqref{fact:extension} to $G_x$ and the Zariski closure
  of $\mathscr{W}$ shows that
  $\pi(\mathscr{W})$ is also contained in the set of points which satisfy {\em (f)}.
  To prove the reverse inclusion, it suffices to apply
  Fact~\ref{fact:elim}\eqref{fact:extension}
  by lifting the solutions of {\em (a)} from the $\uline{z}$-space to points in
  $W$.

  It remains to show that the extra conditions {\em (b)}, {\em (d)} and {\em
    (e)} then do satisfy $\cfiber{i}{u}{W}$. Note that condition {\em (d)} is a
  subset of {\em (e)} and that the combination of {\em (b)} and {\em (d)} is
  equivalent to saying that for all ${{\bm{\alpha}}}$, all polynomials of degree at most $i-1$ 
  are identically zero (when $\uline{z}$ is specialized to
  ${{\bm{\alpha}}}$) while, thanks to condition {\em (d)} and
  \cite[Chap. 3, \S 5, Theorem 2, p. 156]{CoLiOS07}, we
  can conclude that ${{\bm{\alpha}}}$ has a preimage by $\pi$ in $\mathscr{W}$ and the gcd of
  all other polynomials of higher degree in $i$ in $G_x$ (still with $\uline{z}$
  instantiated to ${{\bm{\alpha}}}$) has degree at least $i$ and has at least $i$
  distinct roots.
\end{proof}

Hence, the algorithm which relies on \Cref{prop:elim} consists in:
\begin{enumerate}
  \item
  Computing a Gröbner basis~$G$ for $\mathcal{J}\subset \field{t}[x, u, \uline{z}]$
  w.r.t. some lexicographical ordering $\succ$ with $x \succ u \succ
  \underline{z}$;
\item Computing the relations that define $\mathcal{S}^{(i)}$;
  \item
  For each conjunction of constraints defining $\mathcal{S}^{(i)}$: eliminating~$u$
  from $G_x$ and from the defining equations
  and~eliminating the saturation variables introduced to
  handle inequations (still, the inequations should be kept in the output).\\
  \end{enumerate}\vspace{-0.2cm}
  Note that in step ($3$), Gröbner bases can be used to perform the elimination.
  This also has the advantage to determine if there are points which do satisfy
  both the equations and inequations defining~$\mathcal{S}^{(i)}$, hence deciding its
  emptiness.\\
  It should be noted that, in practice, one can avoid to use
  $\field{t}$ as a base field and perform the computations in $\mathbb{K}[x, u,
  \uline{z}, t]$ with an elimination ordering where $t$ is the smallest
  variable. Specialization properties of Gröbner bases \cite{Kalkbrener97} show
  that one obtains this way a nonreduced Gröbner basis for $\mathcal{J}$.
  Computing the conditions {\em (a)}-{\em(e)} is straightforward with standard
  computer algebra orderings \cite{GaGe13}.

\smallskip 
For the application to DDEs, we now make the following hypothesis which will be assumed to hold in~\cref{sec:elimination,sec:stickelberger}:
      %
\begin{align*}
\label{hyp3} \tag{\bf H3}
\underline{\textsf{Hypothesis 3:}}\;\;\;\; \cfiber{k}{u}{\overline{\mathcal{X}}}
\text{ and } \mathcal{S}_k(\check{z_1}, \overline{\mathcal{X}}) \text{ are finite sets.}
\end{align*}
      
\begin{proposition}\label{prop:elim_DDE}
Let~$P$ be as in~\eqref{eqn:initpoleqn} and~$a\in\mathbb{K}$.
Assume~\eqref{hyp1}, \eqref{hyp3} and that $P$ is squarefree.
We set~$\mathcal{J} :=\langle P, \partial_x P, \partial_uP\rangle:(u-a)^\infty\subset
\mathbb K(t)[x, u, \underline{z}]$.
If $\mathcal{G}$ is the set of polynomials computed by the algorithm
underlying~\Cref{prop:elim}, then any element $R\in\langle
\mathcal{G}\rangle \cap\mathbb{K}[t, z_0]$
satisfies~$R(t, F(t, a))~=~0$.
\end{proposition}
\begin{proof}
  First, it results from Sard's lemma~\cite[Prop.~$B.2$]{SaSc17}
  that if~$P$ is squarefree, then
  $\mathcal{J}\cap\mathbb{K}(t)[\underline{z}]$ is not reduced to~$0$
  (by using the same proof as in~\cite[Lemma~$2.3$]{BoChNoSa22}, for~$k>1$).
  Now as~$\mathcal{G}$
  is the output (for~$i=k$) of the algorithm on which~\cref{prop:elim} relies, we have
  by~\cref{prop:elim} that the zero set of~$\mathcal{G}$ is the Zariski closure
  of~$\mathcal{F}_k(u, \overline{\mathcal{X}})$. Finally, using~\eqref{hyp1}
  implies that the projection of~$\mathcal{F}_k(u,  \overline{\mathcal{X}})$ onto the
  $z_0$-coordinate space contains the value~$F(t, a)$. Hence using Fact~\ref{fact:elim}\eqref{fact:elimination},
  any element $R\in\langle \mathcal{G}\rangle
  \cap\mathbb{K}[t, z_0]$ satisfies~$R(t, F(t, a))=0$.
\end{proof}

\section{Geometric approach}\label{sec:stickelberger}
\def\mathscr{W}{\mathscr{W}}
\def\mathsf{denom}{\mathsf{denom}}


Let~$\mathscr{V}\subset\acfield{t}^{k+2}$ be an algebraic set
associated to an
ideal~$\mathcal{J}$ of
$\mathbb K(t)[x, u, \underline{z}]$.
For a set of variables (or scalars)
$z_0, \ldots, z_{k-1}$, we define the notation $\check{z_1} := z_0, z_2, \ldots, z_{k-1}$. Also,
we consider the canonical inclusion~$j:\mathbb{K}(t)[x, u, \underline{z}]\rightarrow
\mathbb K(t, \check{z_1})[x, u, z_1]$, where the choice $z_1$ is arbitrary.

In this section, we say that assumption~\hypertarget{hypS}{$\bf{(S)}$} holds
if the following assumptions hold:
\begin{itemize}\label{hypS}
  \item \eqref{hyp:stick} holds,
\item  the restriction~$\pi_{\check{z_1}}|_{\mathscr{V}}$ is
  \emph{dominant} (i.e. $\overline{\pi_{\check{z_1}}(\mathscr{V})}
        =
        \acfield{t}^{k-1}$),
\item $\mathcal{J}$ has dimension~$k-1$ in~$\mathbb K(t)[x,u, \underline{z}]$,
\item $\mathcal{J}_{z_1}:=\langle j(\mathcal{J}) \rangle$ has dimension~$0$
  in~$\mathbb K(t, \check{z_1})[x, u, z_1]$.
\end{itemize}
\indent \;\;\;
Recall that using~\cref{lemma:SW_constructible}, the set~$\mathcal{S}_k(\check{z_1}, \mathscr{V})$ is constructible.
 In this section, we design an algorithm with the following
specification: it takes as input a finite set of polynomials
of~$\mathbb K(t)[x, u, \underline{z}]$ generating a radical ideal $\mathcal{J}$
satisfying assumption~~\hyperlink{hypS}{$\bf{(S)}$} and such that
$\mathcal{J}\cap\mathbb K(t)[\underline{z}]$ is principal;
it returns a
finite set of polynomials such that their
vanishing set is the Zariski
closure of~$\mathcal{S}_k(\check{z_1}, \mathscr{V})$, for~$\mathscr{V}\subset\acfield{t}^{k+2}$ the zero set of~$\mathcal{J}$.

\indent We first want to determine algebraic relations which induce a zero set
containing~$\mathcal{S}_k(\check{z_1}, \mathscr{V})$. The ideal~$\mathcal{J}_{z_1}$ having
dimension~$0$, the quotient ring $\mathbb{K}(t, \check{z_1})[x, u,
z_1]/\mathcal{J}_{z_1}$ defines a $\mathbb{K}(t, \check{z_1})$-vector space of
finite dimension. We introduce the multiplication map $m_{z_1}: f\mapsto f\cdot
z_1$ which maps $\mathbb{K}(t, \check{z_1})[x, u, z_1]/ \mathcal{J}_{z_1} $ to
itself, and consider its characteristic polynomial $\xi\equiv\xi_{z_1} \in\mathbb{K}(t,
\check{z_1})[z_1]$. We denote by~$\chi\equiv\chi_{z_1}
\in\mathbb{K}(t)[\underline{z}]$ the numerator
of~$\xi$ and by $\mathsf{denom}(\xi)\in\mathbb{K}(t)[\underline{z}]$ its denominator.
As~$\mathsf{denom}(\xi)$ depends only on~$k~-~1$ variables,~$V(\mathsf{denom}(\xi))$
is a priori defined
in~$\acfield{t}^{k-1}$. Seeing~$\mathsf{denom}(\xi)$ as a polynomial in~$\mathbb K(t)[x, u, \underline{z}]$
(resp.~$\mathbb K(t)[u, \underline{z}]$,~$\mathbb K(t)[\underline{z}]$), its zero set
is~$\acfield{t}^3\times V(\mathsf{denom}(\xi))$ (resp. $\acfield{t}^2\times V(\mathsf{denom}(\xi))$,
$\acfield{t}\times V(\mathsf{denom}(\xi))$).
With a slight abuse of notation, we denote all these sets by~$V(\mathsf{denom}(\xi))$:
the precise definition domain will be
implicitly dependent on the set with which we intersect/take the complement,
etc.\\
\indent Let~$Z :=
(\overline{\mathbb{K}(t)}^{k+2}\setminus V(\mathsf{denom}(\xi))) \cap
\overline{\mathbb{K}(t)}^3\times~\pi_{\check{z_1}}(\mathscr{V})\subset\acfield{t}^{k+2}$.

\begin{lemma}\label{lemma:correspondance_multmap}
  Under~~\hyperlink{hypS}{$\bf{(S)}$} and assuming the radicality of~$\mathcal{J}$,
  the set~$Z$ is a dense subset of~$\acfield{t}^{k+2}$ which satisfies
$V(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}] +
\langle \chi, \partial_{z_1} \chi, \ldots, \partial_{z_1}^{k-1}\chi\rangle)\cap{\pi(Z)}
= {\mathcal{S}_k(\check{z_1}, \mathscr{V}\cap Z)}.$
\end{lemma}
\begin{proof}
  The density of~$Z$ follows from 2 facts: (i) the complement of~$V(\mathsf{denom}(\xi))$
  is dense in~$\acfield{t}^{k-1}$; (ii) by~~\hyperlink{hypS}{$\bf{(S)}$},
  $\pi_{\check{z_1}}|_{\mathscr{V}}$ is
  dominant.\\
  \indent Now, from the~$0$-dimensionality of~$\mathcal{J}_{z_1}$, the
  eigenvalues of the endomorphism $m_{z_1}$ of
  $\mathbb{K}(t, \check{z_1})[x, u, z_1]/\mathcal{J}_{z_1}$
  are, by~Fact~\ref{fact:elim}\eqref{fact:eigenvalue},
  the $z_1$-coordinates of the zero set~$\mathscr{V}_{z_1}\subset\overline{\mathbb K({t, \check{z_1}})}^3$
  associated to~$\mathcal{J}_{z_1}$. Also note that we can specialize
  $\underline{z}$ in~$\xi$ to any point of $\pi(\mathscr{V}\cap
  Z)$.\\
  \indent We now prove the direct inclusion. We pick ${{\bm{\alpha}}}=(\alpha_0,
  \ldots, \alpha_{k-1})\in V(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}] +
  \langle \chi, \partial_{z_1} \chi, \ldots,
  \partial_{z_1}^{k-1}\chi\rangle)\cap\pi(Z)$, and we consider the specialized
  ideal~$\mathcal{J}_{{\bm{\alpha}}}$ obtained by specialization of
  $\mathcal{J}$ to~$\check{z_1}=\check{{{\bm{\alpha}}}}$. As
  $(\check{{{\bm{\alpha}}}})\in
  \pi_{\check{z_1}}(Z)\subset\pi_{\check{z_1}}(\mathscr{V})$, it follows that~$\mathcal{J}_{{\bm{\alpha}}}$
  has dimension~$0$ and that its zero
  set corresponds to the zero set of the specialized ideal~$\mathcal{J}_{z_1, {{\bm{\alpha}}}}$
  obtained after specialization of~$\mathcal{J}_{z_1}$ to
  $\check{z_1}=\check{{{\bm{\alpha}}}}$. As
  ${{\bm{\alpha}}}\in V(\chi, \partial_{z_1} \chi, \ldots, \partial_{z_1}^{k-1}\chi)$, we
  have that $z_1 = \alpha_1$ is a root of multiplicity
  at least~$k$ of the polynomial $\chi(z_1, \check{z_1} = \check{{{\bm{\alpha}}}})$.
 By the
 radicality of~$\mathcal{J}$, and the use of~Fact~\ref{fact:elim}\eqref{fact:eigenvalue}
 in the
  $0$-dimensional ideal $\mathcal{J}_{z_1, {{\bm{\alpha}}}}$, there exist by the correspondence
  $V(\mathcal{J}_{{\bm{\alpha}}}) = V(\mathcal{J}_{z_1, {{\bm{\alpha}}}})$ at least
  $k$ points in~$\mathscr{V}\cap(\cap_{i=0, i \neq 1}^{k-1} \{z_i=\alpha_i\})$.
  Hence~$\alpha\in\mathcal{S}_k(\check{z_1}, \mathscr{V}\cap Z)$.
    
  \indent We now prove the reverse inclusion. Let~${{\bm{\alpha}}}=(\alpha_0, \ldots,
  \alpha_{k-1}) \in\mathcal{S}_{k}(\check{z_1}, \mathscr{V} \cap Z)$. By definition
  of~$\mathcal{S}_{k}(\check{z_1}, \mathscr{V}\cap Z)$, we have that ${{\bm{\alpha}}}\in\pi(\mathscr{V}\cap Z)$.
  As we have~$\pi(\mathscr{V}\cap Z)\subset\pi(\mathscr{V})\cap\pi(Z)$, it follows that
  ${{\bm{\alpha}}}\in\pi(Z)$. It remains to show that~${{\bm{\alpha}}}\in
  V(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}] + \langle \chi, \ldots,
  \partial_{z_1}^{k-1}\chi\rangle)$. By Fact~\ref{fact:elim}\eqref{fact:elimination},
  we have that~${{\bm{\alpha}}}\in
  V(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}])$. Now,
  as~$\pi_{\check{z_1}}({{\bm{\alpha}}})\notin V(\mathsf{denom}(\xi))$, we can consider for $0\leq i
  \leq k-1$ the well-defined polynomials $\partial_{z_1}\chi (z_1, \check{z_1} =
  \check{{{\bm{\alpha}}}})$. By the set equality $V(\mathcal{J}_{{\bm{\alpha}}}) = V(\mathcal{J}_{z_1, {{\bm{\alpha}}}})$,
  the fact that there exist
  at least $k$ points in~$\mathscr{V}\cap(\cap_{i=0, i\neq1}^{k-1}\{z_i=\alpha_i\})$
  translate, by the application of Fact~\ref{fact:elim}\eqref{fact:eigenvalue}
  to the
  well-defined $0$-dimensional ideal~$\mathcal{J}_{z_1, {{\bm{\alpha}}}}$,
  to the vanishings
  ~$\partial_{z_1}^i\chi(\underline{z}={{\bm{\alpha}}})=0$, for all~$0\leq i \leq k-1$.
\end{proof}

In all the combinatorial examples that we considered so far,
taking $\mathcal{J} := \langle P, \partial_x P, \partial_u P\rangle:(u-a)^\infty$
for~$P$ as in~\eqref{eqn:initpoleqn}
always led to~$\mathcal{S}_k(\check{z_1}, \overline{\mathcal{X}}\cap Z) =
\mathcal{S}_k(\check{z_1}, \overline{\mathcal{X}})$.
 %
Hence in order to simplify things and to avoid introducing additional technicalities, 
we assume in the rest of this section that
\begin{align*}
  \label{hyp:ua} \tag{${\mathcal{Z}}$}
  \text{ $\mathcal{S}_k(\check{z_1}, \mathscr{V}\cap Z) = \mathcal{S}_k(\check{z_1}, \mathscr{V})$.}
\end{align*}
From an application viewpoint, working under this new generic assumption is
(as mentioned above) harmless. 

  \begin{lemma}\label{lemma:characterization_stickelberger}
    Assuming~\eqref{hyp:ua},~\hyperlink{hypS}{$\bf{(S)}$},
    and the radicality of~$\mathcal{J}$,
    the Zariski closure of~$\mathcal{S}_k(\check{z_1}, \mathscr{V})$ is the zero set
    of the saturated ideal
    $(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}]+
    \langle \chi, \ldots, \partial_{z_1}^{k-1}\chi\rangle):d(\xi)^\infty$.
  \end{lemma}
  \begin{proof}
    We define $V_1:=V(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}]+
    \langle \chi, \ldots, \partial_{z_1}^{k-1}\chi\rangle)$.
    Using~\eqref{hyp:ua} we can
    replace~$\mathcal{S}_k(\check{z_1}, \mathscr{V})$ by~$\mathcal{S}_k(\check{z_1}, \mathscr{V}\cap Z)$.
    Also, \cref{lemma:correspondance_multmap} implies that
    the Zariski closure of~$\mathcal{S}_k(\check{z_1}, \mathscr{V})$ is equal to the Zariski closure
    of~$V_1\cap\pi(Z)$. Denote $W:=\overline{\mathbb{K}(t)}\times\pi_{\check{z_1}}(\mathscr{V})\subset
    \overline{\mathbb{K}(t)}\times\overline{\mathbb{K}(t)}^{k-1}$
    (which is dense
    in~$\overline{\mathbb{K}(t)}^{k}$ as~$\pi_{\check{z_1}}|_{\mathscr{V}}$ is assumed dominant).
    By the definition of~$Z$, we
    have~$\overline{V_1\cap\pi(Z)} =
    \overline{V_1\cap\big(W\setminus \pi(V(d(\xi)))\big)}$. Using the density
    property of~$W$, we thus have that $\overline{V_1\cap\pi(Z)} =
    \overline{V_1\setminus \pi(V(d(\xi)))}$. The results hence follows
    as the zero set of~$(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}]+
    \langle \chi, \ldots, \partial_{z_1}^{k-1}\chi\rangle):d(\xi)^\infty$ is
    $\overline{V_1\setminus \pi(V(d(\xi)))}$.
\end{proof}


Applying~\cref{lemma:correspondance_multmap,lemma:characterization_stickelberger}
and under~\eqref{hyp:ua}, our aim is hence
to compute~$\mathcal{J}\cap\mathbb{K}(t)[\underline{z}] +
\langle \chi, \partial_{z_1} \chi, \ldots, \partial_{z_1}^{k-1}\chi\rangle$.
Regarding complexities estimates,
the use of parametric geometric resolution tools (\cite{Schost03}) allows to prove
the following result.
\begin{lemma}\label{lemma:complexity_stickelberger}
  Assume that~$\mathcal{J}$ is radical
  of degree~$D$, of dimension~$k-1$ and
  that: $\mathcal{J}\cap\mathbb{K}(t)[\check{z_1}]$ is principal and
  $\mathcal{J}$
  induces an ideal~$\mathcal{J}_{z_1}\subset\mathbb{K}(t, \underline{z})[x,u, z_1]$
  of dimension~$0$.
Suppose given a straight-line program of length~$L$ evaluating the polynomials
defining~$\mathcal{J}$. 
Then computing a generator
of~$\mathcal{J}\cap\mathbb{K}[t, \underline{z}]$
and the polynomials $\chi, \partial_{z_1}\chi, \ldots,
\partial_{z_1}^{k-1}\chi$ is done
in~$\tilde{O}(D^{k+1}\cdot 8^k\cdot(L+k^2) + k\cdot D^{2(k+1)})$ ops. in~$\mathbb{K}$.
\end{lemma}
\begin{proof}
We first denote by~$\mathcal{J}_{z_1}$ the induced ideal of dimension~$0$
in~$\mathbb{K}(t, \check{z_1})[x,u, z_1]$.
Using the radicality and the dimension assumptions, it is possible
to compute a parametric geometric resolution of the zero set of~$\mathcal{J}_{z_1}$.
Applying the algorithm underlying~\cite[Theorem~$2$]{Schost03},
we compute $A, B\in\mathbb{K}(t, \hat{\underline{z}})[\lambda]$
which give a parametrization $z_1= A(\lambda)/\partial_\lambda B(\lambda)$ over the field
extension defined by $B(\lambda)=0$. Using~\cite[Theorem~$2$]{Schost03}, this can be done
in~$\tilde{O}(D^{k+1}\cdot 8^k\cdot (L+k^2))$ ops.
in~$\mathbb{K}$. Applying Fact~\ref{fact:elim}\eqref{fact:Stickelberger},
$\chi$ is equal to
the resultant in~$\lambda$ of the numerators of $B$
and $z_1\cdot\partial_\lambda B-A$. Using now fast
computation of bivariate resultants~\cite[\S5]{HyNeSc19} and using the
upper bound~$D$~\cite[Theorem~$1$]{Schost03} on the partial degrees of the coefficients
of both the numerator and denominator of $A$ and~$B$, we obtain that the number of
evaluation-interpolation points needed for $t, \check{z_1}$ is
$O(D^{2k})$. This gives a cost in~$\tilde{O}(D^{2k+1.63})$ ops. in~$\mathbb{K}$ for computing
the numerator of $\chi$, whose partial degrees are bounded by~$D^2$. Finally,
if $L'$ is the length of a straight-line program evaluating~$\chi$, 
Theorem~1 in \cite{BaSt83} allows us to compute
$\partial_{z_1}\chi, \ldots, \partial_{z_1}^{k-1}\chi$ using $\tilde{O}(k\cdot L')
\subset \tilde{O}(k\cdot D^{2(k+1)})$
ops. in~$\mathbb{K}$. This yields a final complexity
in~$\tilde{O}(D^{k+1}\cdot 8^k\cdot(L+k^2) + k\cdot D^{2(k+1)})$ ops. in~$\mathbb{K}$.
As~$\mathcal{J}$ is assumed radical,
a generator of~$\mathcal{J}\cap\mathbb{K}[t, \underline{z}]$ is given
by the squarefree part of~$\chi$. 
The cost of this squarefree part computation is negligible and absorbed in
the above complexity.
\end{proof}


\begin{remark}\label{rmk:key_geom_sat}
{(${\mathbf i}$)} When the characteristic polynomial of $m_{z_1}$ equals its minimal polynomial,
 \cref{lemma:complexity_stickelberger}
amounts to computing a generator of~$\mathcal{J}\cap\mathbb{K}[t, \underline{z}]$.
In this case, the complexity of \cref{lemma:complexity_stickelberger} drops
to~$\tilde{O}(D^{k+1}\cdot 8^k\cdot(L+k^2) + k\cdot D^{(k+1)})\subset
\tilde{O}(D^{k+1}\cdot 8^k\cdot(L+k^2))$. \\
{(${\mathbf {ii}}$)}
\cref{lemma:complexity_stickelberger}
allows with the same complexity to compute the characteristic polynomial~$\xi_u$ of
the multiplication map $m_u$. Denoting by~$\chi_u$ the numerator of~$\xi_u$,
the refinement of $\mathcal{S}_k(\check{z_1}, \mathscr{V})$ consisting in
counting only the distinct solutions w.r.t.~$u$
is equivalent to considering the polynomial
$m\cdot \textsf{disc}_u(\chi_u)-1$ (and then eliminating~$m$).
Another useful practical refinement consists in adding to the polynomials
defining~$\overline{\mathcal{S}_k(\check{z_1}, \mathscr{V})}$
all the polynomials defining~$\overline{\mathcal{S}_k(\check{z_i}, \mathscr{V})}$,
for all~$0\leq i \leq k-1$ (by a slight adaptation of the present section).
\end{remark}


\cref{lemma:complexity_stickelberger,lemma:characterization_stickelberger}
yield an algorithm and prove its complexity.
As most of the combinatorial examples we have encountered
so far are stated with~$\mathbb{K}=\mathbb{Q}$, our aim is to 
 take in practice the benefit of fast multi-modular arithmetic. Also, we prefer to reduce
 the computations in~$\mathbb{Q}$ instead of~$\mathbb{Q}(t)$ by using evaluation-interpolation
 on the parameter~$t$.
 In practice, the underlying algorithm
 makes use of the specialization properties of
Gr\"obner bases~\cite[Prop.~$1$, p.~$308$]{CoLiOS07} (see~\cite[\S4.4]{LeSa22} for details).


\begin{proposition}\label{prop:stick_DDE}
Let~$P\in\mathbb{K}[x, \underline{z}, t, u]$
be as in~\eqref{eqn:initpoleqn}, $a\in\mathbb{K}$ and
define $\mathcal{J}:=\langle P, \partial_xP, \partial_uP\rangle :(u-a)^\infty
\subset\mathbb{K}(t)[x,u, \underline{z}]$.
We assume that:~\eqref{hyp1}, \eqref{hyp3}, ~\eqref{hyp:ua} and~\hyperlink{hypS}{$\bf{(S)}$}
hold and that
the input assumptions of~\cref{lemma:complexity_stickelberger}
are satisfied by the ideal
$\mathcal{J}$.
Then
any nonzero $R\in(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}] +
\langle \chi, \partial_{z_1} \chi, \ldots, \partial_{z_1}^{k-1}\chi\rangle)\cap
\mathbb{K}[t, z_0]$ satisfies~$R(t, F(t, a))=0$.
\end{proposition}
\begin{proof}
  Using the algorithm underlying the proof of~\cref{lemma:complexity_stickelberger},
  we compute a generator of
  $\mathcal{J}\cap\mathbb{K}(t)[\underline{z}]$ and 
  $\chi, \partial_{z_1} \chi, \ldots, \partial_{z_1}^{k-1}\chi$.
  Now using~\cref{lemma:characterization_stickelberger},
  the zero set of~$(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}] +
  \langle \chi, \partial_{z_1} \chi, \ldots, \partial_{z_1}^{k-1}\chi\rangle)$
  is~$\mathcal{S}_k(\check{z_1}, \overline{\mathcal{X}})$.
  Hence applying Fact~\ref{fact:elim}\eqref{fact:elimination},
  the zero set of $(\mathcal{J}\cap\mathbb{K}(t)[\underline{z}] +
\langle \chi, \partial_{z_1} \chi, \ldots, \partial_{z_1}^{k-1}\chi\rangle)\cap
\mathbb{K}[t, z_0]$ is the Zariski closure of the
projection of $\mathcal{S}_k(\check{z_1}, \overline{\mathcal{X}})$ onto the $z_0$-coordinate space, which by~\eqref{hyp1} contains~$F(t, a)$. Hence
as the latter elimination ideal is not reduced to 0
(since~$\mathcal{S}_k(\check{z_1}, \overline{\mathcal{X}})$
is finite), any element of it annihilates~$F(t, a)$.
\end{proof}



\section{Conclusion and perspectives}\label{sec:conclusion}

We have introduced new fast algorithms for solving arbitrary order DDEs with one catalytic variable.
These algorithms are not only fast in theory, but also efficient in practice; moreover, they allow to
solve nontrivial combinatorial applications, as showed by the experimental results in
\cref{sec:experiments}. 
Our implementations yielded practical improvements for a large majority of them, and allowed us to solve
one (5-constellations)
on which the state-of-the-art methods could not terminate. 

For future works, we wish to develop and to provide complete implementations of the
algorithms that we have introduced in the present paper, and to make them available to the combinatorics
community.
A different, more theoretical direction, is to
pursue the geometrical investigation analysis of the problem of computing~\emph{exceptional fibers}
initiated in~\cref{sec:elimination,sec:stickelberger}.



\section{Experiments}\label{sec:experiments}

  \paragraph{Aim} We first report on practical variants of the hybrid
  guess-and-prove method \textsf{(hgp)} and then provide and analyze tables
  of our implementations of
  \cref{sec:direct,sec:guessandprove,sec:elimination,sec:stickelberger}.
  %
   The benchmark DDEs on which we test our various implementations 
  have combinatorial origins and the literature qualifies their
  resolution as a highly non-trivial problem. 
  More precisely, we consider solving:
  Eq.~(4.22) in \cite{Bernardi08}
  (near-triangulations), 
  Prop.~12 in \cite{BMJ06} ($m$-constellations, $m\in\{4, 5\}$) 
  and
  Eq.~(3) in \cite{BMFP12} ($3$-Tamari lattices).\\
  \indent Recall that the input of the algorithms
  in~\cref{sec:direct,sec:guessandprove,sec:elimination,sec:stickelberger} 
  consists of 
  a polynomial
  $P$ as in~\eqref{eqn:initpoleqn} and 
  of
  a specialization point $a\in\mathbb{K}$, while
  their output
  is, under~\eqref{hyp3} and up to eliminating variables,
  a non-zero polynomial~$R\in\mathbb{K}[t, z_0]$
  such that~$R(t, F(t, a)) = 0$.
  %
  %
  \paragraph{Implementations}
  The DDEs we study being defined over~$\mathbb K=\mathbb Q$, we use multi-modular arithmetic
  and CRT (Chinese remainder theorem) for
  \cref{sec:direct,sec:elimination,sec:stickelberger}.
  Also, we reduce the computation from~$\mathbb Q(t)$ to
  $\mathbb Q$ by performing ev.-int. (evaluation-interpolation) on either~$t$ or~$z_0$.
  We incorporate~ $\mathbf{(ii)}$ of Rmk. \ref{rmk:key_geom_sat}, but do
  not use the inequalities describing~$\mathcal{S}^{(k)}$ in~\cref{sec:elimination}. 
  Finally,
  we use standard tools in computer algebra to improve each of our
  implementations.\\
  \indent The practical variant
  of the (\textsf{hgp}) strategy mentioned above has the following
  motivation.
When one performs (say, for a ``random'' prime~$p$) the computation of a modular image~$R_p$
of~$R$ by using one of~\cref{sec:direct,sec:elimination,sec:stickelberger},
the computation (if it ends) gives access to the partial degrees
of~$R\in\mathbb{Q}[t, z_0]$. If either one modular computation or
the lift over~$\mathbb{Q}$ is too time consuming,
the following variant of \textsf{hgp} exploits the knowledge of those partial degrees:
\begin{itemize}
\item Pick a random prime~$2^{27} < p < 2^{31}$ and
   $\theta\in\{1, \ldots, p-1\}$,
\item
  Compute $R_p(t, \theta) \in\mathbb{F}_p[t]$ and $R_p(\theta, z_0)\in
  \mathbb{F}_p[z_0]$ using one
  of~\cref{sec:direct,sec:elimination,sec:stickelberger};
  set $d_t$ and $d_{z_0}$ their respective degrees,
\item Compute~$F(t, a)\mod t^{2d_td_{z_0}+1}$,
\item Guess $\tilde{R}\in\mathbb{Q}[t, z_0]$ s.t.
  $\tilde{R}(t, F(t, a))=O(t^{(d_t+1)(d_{z_0}+1)-1})$,
\item Check that~$\tilde{R}(t, F(t, a)) = O(t^{d_t\cdot \deg_{z_0}(\tilde{R}) +
  \deg_t(\tilde{R})\cdot d_{z_0}+1})$.
\end{itemize}
The above algorithm is a simple extension of the one in~\cite[\S2.2.2]{BoChNoSa22},
with the total degree replaced by partial degrees. Also, our implementation generates
terms of~$F(t,a)$ by first computing terms of~$F(t, u)$ and then specializing to~$u=a$.
Any optimization of this step would result in much better timings for
the~$\textsf{(hgp)}$.


In our experiments, we consider the following data:
\begin{itemize}
  \item $\mathbf{S}$: sections (and hence algorithms) used,
\item $\mathbf{\#\mathcal{P}}$: number of primes used for the CRT,
  \item $\mathbf{Z}\in\{t, z_0\}$ the variable on which we perform ev.-int.,
  \item \textbf{$\mathbf{\#}$pts}: number of ev.-int. points needed in~$Z$,
  \item $\mathbf{d}_{\textbf{cp}}$: critical pairs of maximal degree in GB computations,
    \item $\mathbf{d}_\mathbf{M}$: Macaulay matrix maximal size in GB computations (F$4$),
      \item $\mathbf{d}_{\mathcal{I}_\infty}$: degree of the ideal~$\mathcal{I}_\infty$
    in~\cref{sec:direct},
  \item $\mathbf{d}_\chi$: degree of~$\chi_{z_1}\in\mathbb{Q}[t, \underline{z}]$,
  \item $\mathbf{T}$: total timing needed to obtain an output in~$\mathbb{Q}[t, z_0]$,
  \item $\mathbf{d}_\mathbf{Z}$: degree in~$Z$ of output~$R\in\mathbb{Q}[t, z_0]$ s.t.
    $R(t, F(t, a))=0$,
\item $\mathbf{\sigma}$: truncation order in the expansion of~$F(t, a)$,
\item $\mathbf{G}$:
  time spent for guessing an annihilating polynomial in~$\mathbb{Q}[t, z_0]$,
\item $\mathbf{P}$: time spent to prove the guess.
\end{itemize}
\indent \;\;\;The timings are given in seconds (s.), minutes (m.), hours (h.) and days (d.).
The symbol $\infty$ (resp. $-$) means
that the computation (resp. the data) did not finish (resp. was not known)
after 5 days.\\% of computation.\\
\indent All computations
were conducted using Maple on a computer equipped with
Intel® Xeon® Gold CPU 6246R v4 @ 3.40GHz and 1.5TB of RAM with a single thread.
All Gr\"obner bases computations were performed using the  \textsf{C} library
\href{https://msolve.lip6.fr/}{msolve}~\cite{msolve}, 
and all guessing computations were performed using the
\href{http://perso.ens-lyon.fr/bruno.salvy/software/the-gfun-package/}{gfun} Maple
package~\cite{Gfun}.

We obtain the following tables:

\tiny
\begin{center} {\bf{\footnotesize \cite[Proposition~$\mathbf{4.3}$]{Bernardi08}}, {\bf $k=\mathbf{2}$}}\\
  \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}
    \hline
    $\mathbf{S}$ & $\mathbf{\#\mathcal{P}}$ & $\mathbf{Z}$ & \textbf{$\mathbf{\#}$pts}
    & \textbf{$\mathbf{d}_{\textbf{cp}}$} & \textbf{$\mathbf{d}_\mathbf{M}$} &
    \textbf{$\mathbf{d}_{\mathbf{\mathcal{I}}_\infty}$} &
    \textbf{$\mathbf{d}_{\mathbf{\chi}}$}
    & \textcolor{blue}{$\mathbf{T}$} & $\mathbf{d}_\mathbf{t}$ & $\mathbf{d}_{\mathbf{z_0}}$\\
    \hline
    \cref{sec:direct} & $3$ & $t$ & $265$ & $18$ & $4\cdot10^5\times4\cdot10^5$ & $12$ & $\times$ & \textcolor{blue}{$55$s} & $132$ & $6$\\
    \cref{sec:elimination}  & $3$ & $t$ & $265$ & $206$ & $2\cdot10^4\times2\cdot10^4$ &$\times$  & $\times$ & \textcolor{blue}{$1$m$10$s} & $132$ & $6$\\
    \cref{sec:stickelberger} & $\times$ & $\times$ & $\times$ & $\times$ &
    $\times$
    & $\times$ & $469$ & \textcolor{blue}{$30$s} & $1173$ & $33$ \\
    \hline
\end{tabular}
\end{center}

\tiny
\begin{center}
  \begin{tabular}{|c|c|c|c|c|c|c|}
    \hline
    $\mathbf{S}$ & $\mathbf{d}_\mathbf{t}$ & $\mathbf{d}_{\mathbf{z_0}}$ &$\mathbf{\sigma}$&
    $\mathbf{G}$  &
  $\mathbf{P}$ & \textcolor{blue}{$\mathbf{T}$} \\
  \hline
  \cref{sec:guessandprove,sec:stickelberger} &  $132$ \textcolor{blue}{($7$m)}
  & $6$ \textcolor{blue}{($0.4$s)}& $\mod t^{2048}$ \textcolor{blue}{($\infty$)} & $-$ & $-$ &
  $-$
 \\

  \hline
\end{tabular}\\[0.2cm]
\end{center}


\tiny
\begin{center}  {\bf {\footnotesize $\mathbf{4}$-constellations}, {\bf $k=\mathbf{3}$}}\\
  \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}
    \hline
    $\mathbf{S}$ & $\mathbf{\#\mathcal{P}}$ & $\mathbf{Z}$ & \textbf{$\mathbf{\#}$pts}
    & \textbf{$\mathbf{d}_{\textbf{cp}}$} & \textbf{$\mathbf{d}_\mathbf{M}$} &
    \textbf{$\mathbf{d}_{\mathbf{\mathcal{I}}_\infty}$} &
    \textbf{$\mathbf{d}_{\mathbf{\chi}}$}
    & \textcolor{blue}{$\mathbf{T}$} & $\mathbf{d}_\mathbf{t}$ & $\mathbf{d}_{\mathbf{z_0}}$\\
    \hline
    \cref{sec:direct} & $3$ & $t$ & $7$ & $9$ & $3\cdot10^4\times3\cdot10^4$ & $42$
    & $\times$ & \textcolor{blue}{$4$m} & $3$ & $7$ \\
    \cref{sec:elimination} & $3$ & $t$ & $7$ & $33$ & $7\cdot 10^3\times 11\cdot 10^3$ & $\times$
    & $\times$ & \textcolor{blue}{$41s$} &  $3$ & $7$ \\
    \cref{sec:stickelberger} & $3$ & $t$ & $7$ &  $51$ &
    $2\cdot 10^5\times 2\cdot 10^5$  & $\times$ & $28$ & \textcolor{blue}{$83s$}
    &  $3$ & $7$  \\
    \hline
\end{tabular}
\end{center}

\tiny
\begin{center} 
  \begin{tabular}{|c|c|c|c|c|c|c|}
    \hline
    $\mathbf{S}$ & $\mathbf{d}_\mathbf{t}$ & $\mathbf{d}_{\mathbf{z_0}}$ &$\mathbf{\sigma}$&
    $\mathbf{G}$  &
  $\mathbf{P}$ & \textcolor{blue}{$\mathbf{T}$} \\
  \hline
 \cref{sec:guessandprove,sec:elimination} & $3$  \textcolor{blue}{($1$s)}  & $7$  \textcolor{blue}{($1.5$s)}& $\mod t^{64}$
 \textcolor{blue}{($18$s)} & \textcolor{blue}{$0.03$s} & \textcolor{blue}{$0.001$s} &
  \textcolor{blue}{$20$s}\\
  \hline
\end{tabular}\\[0.2cm]
\end{center}


\tiny
\begin{center}  {\bf {\footnotesize $\mathbf{3}$-Tamari}, {\bf $k=\mathbf{3}$}}\\
  \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}
    \hline
    $\mathbf{S}$ & $\mathbf{\#\mathcal{P}}$ & $\mathbf{Z}$ & \textbf{$\mathbf{\#}$pts}
    & \textbf{$\mathbf{d}_{\textbf{cp}}$} & \textbf{$\mathbf{d}_\mathbf{M}$} &
    \textbf{$\mathbf{d}_{\mathbf{\mathcal{I}}_\infty}$} &
    \textbf{$\mathbf{d}_{\mathbf{\chi}}$}
    & \textcolor{blue}{$\mathbf{T}$} & $\mathbf{d}_\mathbf{t}$ & $\mathbf{d}_{\mathbf{z_0}}$\\
    \hline
    \cref{sec:direct} & $4$ & $t$ & $11$ & $11$ & $3\cdot10^5\times3\cdot10^5$ & $96$ & $\times$ & \textcolor{blue}{$2$d$2$h} &  $5$ &  $16$\\
    \cref{sec:elimination} & $4$ & $z_0$ & $33$ &$64$ & $10^4\times10^4$& $\times$ & $\times$&   \textcolor{blue}{$2$m} &  $5$ &  $16$ \\
    \cref{sec:stickelberger} & $4$ & $t$ & $11$ & $52$ & $10^4\times10^4$ & $\times$ & $31$ &   \textcolor{blue}{$5$m$42$s}
    & $5$ &  $16$\\
    \hline
\end{tabular}
\end{center}
\tiny
\begin{center}
  \begin{tabular}{|c|c|c|c|c|c|c|}
    \hline
    $\mathbf{S}$ & $\mathbf{d}_\mathbf{t}$ & $\mathbf{d}_{\mathbf{z_0}}$ &$\mathbf{\sigma}$&
    $\mathbf{G}$  &
  $\mathbf{P}$ & \textcolor{blue}{$\mathbf{T}$} \\
  \hline
 \cref{sec:guessandprove,sec:elimination} & $5$  \textcolor{blue}{($0.2$s)}  & $16$  \textcolor{blue}{($5$s)}& $\mod t^{256}$
 \textcolor{blue}{($1$h$40$m)} & \textcolor{blue}{$1$s} & \textcolor{blue}{$0.2$s} &
  \textcolor{blue}{$1$h$40$m}\\
  \hline
\end{tabular}\\[0.2cm]
\end{center}


\tiny
\begin{center} {\bf {\footnotesize $\mathbf{5}$-constellations}, {\bf$k=\mathbf{4}$}}
 
  \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}
    \hline
    $\mathbf{S}$ & $\mathbf{\#\mathcal{P}}$ & $\mathbf{Z}$ & \textbf{$\mathbf{\#}$pts}
    & \textbf{$\mathbf{d}_{\textbf{cp}}$} & \textbf{$\mathbf{d}_\mathbf{M}$} &
    \textbf{$\mathbf{d}_{\mathbf{\mathcal{I}}_\infty}$} &
    \textbf{$\mathbf{d}_{\mathbf{\chi}}$}
    & \textcolor{blue}{$\mathbf{T}$} & $\mathbf{d}_\mathbf{t}$ & $\mathbf{d}_{\mathbf{z_0}}$\\
    \hline
    \cref{sec:direct} & $-$ & $-$ & $-$& $-$& $-$& $-$&$\times$ & \textcolor{blue}{$\infty$} & $-$ & $-$\\
    \cref{sec:elimination}  & $-$ & $t$ & $53$ & $70$ & $2\cdot10^6\times2\cdot10^6$
    & $\times$ & $\times$ &
    \textcolor{blue}{$\infty$} & $26$ & $53$\\
    \cref{sec:stickelberger} & $-$ & $-$&$-$ &$-$&$-$& $\times$ & $-$& \textcolor{blue}{$\infty$} & $-$& $-$\\
    \hline
\end{tabular}
\end{center}


\tiny
\begin{center}
  \begin{tabular}{|c|c|c|c|c|c|c|}
    \hline
    $\mathbf{S}$ & $\mathbf{d}_\mathbf{t}$ & $\mathbf{d}_{\mathbf{z_0}}$ &$\mathbf{\sigma}$&
    $\mathbf{G}$  &
  $\mathbf{P}$ & \textcolor{blue}{$\mathbf{T}$} \\
  \hline
 \cref{sec:guessandprove,sec:elimination} & $26$  \textcolor{blue}{($47$m)}  & $53$  \textcolor{blue}{($45$m)}& $\mod t^{256}$
 \textcolor{blue}{($4$h$35$m)} & \textcolor{blue}{$0.03$s} & \textcolor{blue}{$0.02$s} &
  \textcolor{blue}{$6$h$7$m}\\
  \hline
\end{tabular}\\\vspace{0.2cm}
\end{center}

\normalsize
 \paragraph{Interpretation}
  \indent A first natural observation is that the algorithms
  introduced in 
  \cref{sec:direct,sec:guessandprove,sec:elimination,sec:stickelberger}, 
  as well as their practical variants, are relevant in practice. 
  Moreover, for all the examples, there is always one of the
  new methods which is more efficient in terms of timings than the
  state-of-the-art~(\cref{sec:direct}). According to the tables, there is generally no unique
  method that is always better than the others. On the contrary, the experiments
  show that all the new methods can be useful in practice, depending on the DDEs 
  under study
  (and hence on the properties of the associated zero sets).
  \\
  \indent We now explain the tables related to~$5$-constellations. Note first that
  neither~\cref{sec:direct} nor~\cref{sec:stickelberger} allow to compute any single
  specialization (at~$t$ or~$z_0$ specialized) of a modular image of~$R$.
  Now applying~\cref{sec:elimination}, we manage to compute two specializations
  (first in~$t$, then in~$z_0$) of~$R_p$, for some ``random'' prime $p$. This hence
  gives all the relevant data of the line except $\mathbf{\#\mathcal{P}}$ and~$\mathbf{T}$.
  Those two specializations take respectively~$45$m. and~$47$m. The degrees obtained being
  $d_{z_0}=53$ and~$d_t=26$,
  we would need approximately $53\cdot45$min$=39$h. for each modular computation.
  Estimating the number of such modular computations to be~$5$ (which is very likely
  optimistic), we would hence need at least~$8$ days. Instead of this, we
  use the practical variant of~\cref{sec:guessandprove} mentioned in the current section.
  As the degree of the guessed
  polynomial is low ($\deg_t=2$ and $\deg_{z_0} = 5)$,
  it allows us to compute $256$ terms of the series~$F(t,1)$
  (here $a=1$), and to check the guess with the geometric bounds~$d_t, d_{z_0}$ obtained
  previously.


\balance
\bigskip


