Asymptotic Solutions of Polynomial Equations with Exp-Log Coefficients

TL;DR

This paper introduces an algorithm to compute asymptotic approximations of roots for polynomials with exp-log coefficients, providing explicit expressions and a method to identify real roots, supported by implementation and empirical data.

Contribution

It presents a novel algorithm for asymptotic root approximation of polynomials with exp-log coefficients, including a real root identification method.

Findings

Algorithm successfully computes asymptotic roots with explicit exp-log expressions.

Method effectively distinguishes real roots among approximations.

Empirical data demonstrates the algorithm's practical performance.

Abstract

We present an algorithm for computing asymptotic approximations of roots of polynomials with exp-log function coefficients. The real and imaginary parts of the approximations are given as explicit exp-log expressions. We provide a method for deciding which approximations correspond to real roots. We report on implementation of the algorithm and present empirical data.

Tables

Table 1. Table 1. Example 28 .

$m$	$n$	$2$	$3$	$4$	$5$	$6$	$7$
	Time	$0.171$	$0.235$	$0.292$	$0.276$	$0.389$	$0.448$
$1$	Iter	$4$	$6$	$6$	$7$	$6$	$6$
	LC	$56$	$84$	$124$	$195$	$192$	$323$
	Time	$0.231$	$0.508$	$1.18$	$1.20$	$1.02$	$1.87$
$5$	Iter	$20$	$32$	$44$	$59$	$58$	$67$
	LC	$136$	$229$	$498$	$668$	$588$	$1933$
	Time	$0.376$	$0.796$	$3.53$	$5.25$	$1.64$	$5.18$
$10$	Iter	$40$	$63$	$86$	$124$	$118$	$137$
	LC	$256$	$433$	$896$	$1303$	$1188$	$4008$
	Time	$0.968$	$3.21$	$26.5$	$25.4$	$5.47$	$20.1$
$20$	Iter	$80$	$124$	$170$	$254$	$238$	$277$
	LC	$496$	$837$	$1716$	$2573$	$2388$	$8158$

Equations

$$\mathbb{R}\setminus \{0\}\ni x\to 1/x=x\exp (-\log (x^2))\in \mathbb{R}$$

$${\mathbb{R}}_{+}\ni x\to x^r=\exp (r\log (x))\in \mathbb{R}$$

$$|y-y^{\ast }|\le |\frac{5P(x,y)}{\frac{\partial P}{\partial y}(x,y)}|$$

$$\underset{x\to \infty }{\lim }\frac{\log |f(x)|}{\log |g(x)|}$$

$$\underset{x\to \infty }{\lim }\frac{\log |f(x)|}{\log |g(x)|}=0$$

$$\underset{x\to \infty }{\lim }f(x)g(x{)}^e=0$$

$$\underset{x\to \infty }{\lim }\frac{\log |f(x)|}{\log |g(x)|}=0$$

$$\underset{x\to \infty }{\lim }\frac{\log |f(x)|}{\log |g(x)|}\in \mathbb{R}\setminus \{0\}$$

$$p=a_n{z}^n+\dots +a_0=a_n(z-r_1)\cdot \dots \cdot (z-r_n)\in \mathbb{C}[z]$$

$$q=b_m{z}^m+\dots +b_0\in \mathbb{C}[z]$$

$$q=b_m(z-s_1)\cdot \dots \cdot (z-s_m)$$

$$|q(c)-p(c)|\le \sum _{k=0}^n|b_k-a_k||z{|}^k+\sum _{k=n+1}^m|b_k||z{|}^k\le \delta \sum _{k=0}^m\frac{1}{{ϵ}^k}$$

$$\delta \sum _{k=0}^m\frac{1}{{ϵ}^k}<|a_n|\frac{(1-ϵ){ϵ}^{m+n}}{1-{ϵ}^{m+1}}\frac{1-1/{ϵ}^{m+1}}{1-1/ϵ}=|a_n|{ϵ}^n$$

$$b_i:=\{\begin{array}{ll}{\hat{u}}_i& e_{u,i}<e_{v,i}\\ {\hat{u}}_i+ı{\hat{v}}_i& e_{u,i}=e_{v,i}\\ ı{\hat{v}}_i& e_{u,i}>e_{v,i}\end{array}$$

$$\underset{x\to \infty }{\lim }{\omega }^{-(e_i+d)}(a_i^{\uparrow k}-b_i{\omega }^{e_i})=0$$

$${\omega }^{-(e_i+d)}(a_i^{\uparrow k}-b_i{\omega }^{e_i})={\omega }^{-(e_{u,i}+d)}(u_i^{\uparrow k}-{\hat{u}}_i{\omega }^{e_{u,i}})+ı{\omega }^{-(e_{u,i}+d)}{v}_i^{\uparrow k}$$

$$\underset{x\to \infty }{\lim }{\omega }^{-(e_{u,i}+d)}(u_i^{\uparrow k}-{\hat{u}}_i{\omega }^{e_{u,i}})=$$

$$\underset{x\to \infty }{\lim }{\omega }^{\hat{d}-d}{\omega }^{-(e_{u,i}+\hat{d})}(u_i^{\uparrow k}-{\hat{u}}_i{\omega }^{e_{u,i}})=0$$

$$\underset{x\to \infty }{\lim }{\omega }^{-(e_{u,i}+d)}{v}_i^{\uparrow k}=$$

$$\underset{x\to \infty }{\lim }{\omega }^{(e_{v,i}-e_{u,i})+(\hat{d}-d)}{\omega }^{-(e_{v,i}+\hat{d})}(v_i^{\uparrow k}-{\hat{v}}_i{\omega }^{e_{v,i}})+$$

$${\hat{v}}_i{\omega }^{(e_{v,i}-e_{u,i})-d}=0$$

$$0=\underset{x\to \infty }{\lim }\frac{Q(x,\alpha )}{a_n^{\uparrow k}{\alpha }^n}=\underset{x\to \infty }{\lim }1+\sum _{i=1}^n\frac{a_{n-i}^{\uparrow k}}{a_n^{\uparrow k}{\alpha }^i}$$

$$\underset{x\to \infty }{\lim }|\sum _{i=1}^n\frac{a_{n-i}^{\uparrow k}}{a_n^{\uparrow k}{\alpha }^i}|=0$$

$$\underset{x\to \infty }{\lim }|\sum _{i=1}^n\frac{a_{n-i}^{\uparrow k}}{a_n^{\uparrow k}{\alpha }^i}|=\infty$$

$$c_i={\omega }^{-(e_i+d)}(a_i^{\uparrow k}-b_i{\omega }^{e_i})$$

$$\gamma =\underset{x\to \infty }{\lim }\frac{\log |\beta {\omega }^{\gamma }|}{\log |\omega |}=\gamma +\underset{x\to \infty }{\lim }\frac{\log |\beta |}{\log |\omega |}$$

$$Q(x,\alpha )=\sum _{i\in I}(b_i{\beta }^i{\omega }^{e_i+\gamma i}+c_i{\beta }^i{\omega }^{e_i+d+\gamma i})$$

$$\mu <\nu <\min (\mu +d,\underset{i\in I\setminus J}{\min }{e}_i+\gamma i)$$

$$0={\omega }^{-\nu }Q(x,\alpha )={\omega }^{\mu -\nu }\sum _{i\in J}{b}_i{\beta }^i+\sum _{i\in I\setminus J}{b}_i{\beta }^i{\omega }^{e_i+\gamma i-\nu }+\sum _{i\in I}{c}_i{\beta }^i{\omega }^{e_i+d+\gamma i-\nu }$$

$$\underset{x\to \infty }{\lim }{\omega }^{\mu -\nu }\sum _{i\in J}{b}_i{\beta }^i=0$$

$$Q_{\gamma }(x,z)=R_{\gamma }(x,z)+{\omega }^{\nu -\mu }(\sum _{i\in I\setminus J}{b}_i{\omega }^{e_i+\gamma i-\nu }{z}^i+\sum _{i\in I}{c}_i{\omega }^{e_i+d+\gamma i-\nu }{z}^i)$$

$$|b_s|\frac{(1-ϵ){ϵ}^{n+s}}{1-{ϵ}^{n+1}}>|b_s|{ϵ}^{n+s+1}>{ϵ}^{2n+2}=\delta$$

Full text

1. Introduction

Definition 1**.**

The set of exp-log functions is the smallest set of partial functions $\mathbb{\mathbb{R}\rightarrow R}$ containing $\exp$, $\log$, the identity function and the constant functions, closed under addition, multiplication and composition of functions.

The domain $D(f)$ of an exp-log function $f$ is determined as follows:

1. (1)

the domain of $exp$, the identity function and the constant functions is $\mathbb{R}$ and the domain of $log$ is $\mathbb{R}_{+}$, 2. (2)

$D(f+g)=D(fg)=D(f)\cap D(g)$, 3. (3)

$D(f(g))=g^{-1}(D(f))$.

In particular, $D(f)$ is an open set and $f$ is $C^{\infty}$ in $D(f)$.

Remark 2*.*

The multiplicative inverse function

[TABLE]

and the real exponent power functions

[TABLE]

for $r\in\mathbb{R}$, are exp-log functions.

The domain of an exp-log function consists of a finite number of open, possibly unbounded, intervals and an exp-log function has a finite number of real roots. An algorithm computing domains and isolating intervals for real roots of exp-log functions is given in [11, 12].

We say that a partial function $f:\mathbb{R}\rightarrow\mathbb{C}$ is defined near infinity if $D(f)\supseteq(c,\infty)$ for some $c\in\mathbb{R}$.

Definition 3**.**

A Hardy field [1] is a set of germs at infinity of real-valued functions that is closed under differentiation and forms a field under addition and multiplication.

Theorem 4**.**

[4, 3]** The germs at infinity of exp-log functions defined near infinity form a Hardy field.

Let $P(x,y)=a_{n}(x)y^{n}+\ldots+a_{0}(x)$ where, for $0\leq i\leq n$, $a_{i}(x)=u_{i}(x)+\imath v_{i}(x)$ and $u_{i}$ and $v_{i}$ are exp-log functions defined near infinity ($\imath$ denotes the imaginary unit).

Problem 5**.**

Describe the asymptotic behaviour of roots of $P$ in $y$ as $x$ tends to infinity.

The following theorem [8, 9] shows that the problem is well posed, that is the roots of $P$ in $y$ are $C^{\infty}$ functions in $x$ defined near infinity.

Theorem 6**.**

If $H$ is a Hardy field, then there exists a Hardy field $K\supseteq H$ such that $K[\imath]$ is algebraically closed.

Corollary 7**.**

There exists a Hardy field $K$ such that $P$ has $n$ roots in $y$ (counted with multiplicities) in $K[\imath]$.

Definition 8**.**

Let $f:\mathbb{R}\rightarrow\mathbb{C}$ be a partial function defined near infinity. We say that partial functions $f_{1},\ldots,f_{m}:\mathbb{R}\rightarrow\mathbb{C}$ defined near infinity form an $m$-term asymptotic approximation of $f$ if, for $1\leq i<m$, $\lim_{x\rightarrow\infty}\frac{f_{i+1}(x)}{f_{i}(x)}=0$ and $\lim_{x\rightarrow\infty}\frac{f(x)-\sum_{i=1}^{m}f_{i}(x)}{f_{m}(x)}=0$.

In this paper we present an algorithm which computes asymptotic approximations of roots of $P$ in $y$. The approximations are given as exp-log expressions. The algorithm makes use of the theory of “most rapidly varying” subexpressions developed in [2] to compute limits of exp-log functions. In fact our algorithm applies in the more general case of MrvH fields. The algorithm is based on a Newton polygon technique [6, 14, 13] extended to “series” with arbitrary real exponents.

Algorithms given in [10, 13] solve the problem of finding asymptotic solutions of polynomial equations in more general settings. We chose to extend the algorithm of [2] because we find it simpler to implement and we can give a direct and elementary proof that the computed expressions satisfy our (weaker) requirements.

Example 9**.**

Let $P(x,y)=y^{5}-\exp(x)y-\log(x)$. One-term asymptotic approximations of roots of $P$ in $y$ computed with our algorithm are $r_{1}=-\exp(-x)\log(x)$, $r_{2}=-\exp(-x)^{-1/4}$, $r_{3}=\exp(-x)^{-1/4}$, $r_{4}=-\imath\exp(-x)^{-1/4}$, $r_{5}=\imath\exp(-x)^{-1/4}$. Let us estimate the relative error $\varepsilon_{i}=\lvert\frac{r_{i}-r_{i}^{*}}{r_{i}}\rvert$ of the approximations, where $r_{i}^{*}$ is the exact root closest to $r_{i}$. Using the bound

[TABLE]

on the distance from $y$ to the closest root $y^{*}$ of $P(x,y)$, after simplifications valid for $x>1$, we get $\varepsilon_{1}\leq\frac{5\log(x)^{4}}{\exp(5x)-5\log(x)^{4}}$, and for $i\neq 1$, $\varepsilon_{i}\leq\frac{5}{4}\log(x)\exp(-\frac{5}{4}x)$. Both bounds tend to zero as $x$ tends to infinity and are decreasing for $x>\exp(W_{0}(\frac{4}{5}))=1.63\ldots$, where $W_{0}$ is the principal branch of the Lambert W function. Evaluating the bounds at $x=10$ we get $\log_{10}\varepsilon_{1}\leq-19.5$ and $\log_{10}\varepsilon_{i}\leq-4.96$ for $i\neq 1$. For $x=1000$ we get $\log_{10}\varepsilon_{1}\leq-2167$ and $\log_{10}\varepsilon_{i}\leq-541$ for $i\neq 1$. This shows that we can obtain approximations of roots of $P(1000,y)$ to $500$ digits of precision by evaluating the asymptotic approximations. The evaluation takes $2$ ms. For comparison, direct computation of roots of $P(1000,y)$ to $500$ digits of precision takes $68$ ms. Figure 1.1 shows the asymptotic approximations of the real roots of $P(x,y)$ (dashed curves) and the exact roots (solid green curves).

2. Most rapidly varying subexpressions

In this section we give a very brief summary of terminology and facts necessary for formulating Algorithm 14. We will use the algorithm to compute approximations of coefficients of $P$ in terms of a “most rapidly varying” subexpression present in the coefficients. The algorithm is based on the algorithm MrvLimit described in [2]. For a more detailed introduction and proofs of the stated facts see [2].

Definition 10**.**

The set of exp-log expressions with coefficients in a computable field $C\subseteq\mathbb{R}$ is defined recursively as follows:

1. (1)

elements of $C$ and the variable $x$ are exp-log expressions, 2. (2)

if $f$ and $g$ are exp-log expressions, so are $f+g$, $f\cdot g$ and $\frac{f}{g}$, 3. (3)

if $f$ is an exp-log expression and $c\in C$, then $\exp(f)$, $\log(f)$, and $f^{c}$ are exp-log expressions.

Each exp-log expression represents an exp-log function, however the same function may be represented by many different expressions. In the following, when we refer to the domain, point values, and limits of an exp-log expression, we mean the domain, point values, and limits of the corresponding exp-log function.

Definition 11**.**

Let $E_{\infty}(x)$ be the set of exp-log expressions $f$ such that, for some $c\in\mathbb{R}$, $(c,\infty)\subseteq D(f)$ and either $f=0$ as an expression or $f$ is nonzero on $(c,\infty)$.

Remark 12*.*

Note that we exclude from $E_{\infty}(x)$ expressions that are identically zero in a neighbourhood of infinity, but are not explicitly zero e.g. $f=\exp(\log((x-c)^{2})/2)-x+c$. The algorithm ExpLogRootIsolation of [11] can be used to check whether a given exp-log expression is defined near infinity and to detect expressions that are identically zero in a neighbourhood of infinity and replace them with explicit zeros. ExpLogRootIsolation requires a zero test algorithm for elementary constants. Termination of the currently known zero test algorithm relies on Schanuel’s conjecture [7, 12].

The germs at infinity of functions represented by elements of $E_{\infty}(x)$ form the Hardy field of exp-log functions defined near infinity.

Theorem 13**.**

If $f$ and $g$ are nonzero elements of a Hardy field, then the limit

[TABLE]

exists (in $\overline{\mathbb{R}}=\mathbb{R}\cup\{-\infty,\infty\}$). Moreover, if $\lim_{x\rightarrow\infty}g(x)=0$ and

[TABLE]

then, for any $e>0$,

[TABLE]

The theorem follows from the results in section 3.1.2 of [2].

Following [2], we say that $g$ is more rapidly varying than $f$, or $g$ is in a higher comparability class than $f$, if

[TABLE]

and we denote it $f\prec g$. We say that $f$ and $g$ have the same order of variation, or $f$ and $g$ are in the same comparability class, if

[TABLE]

and we denote it $f\asymp g$. We will also use $f\preceq g$ to denote $f\prec g\vee f\asymp g$. $\omega$ is a most rapidly varying subexpression of $f$ if $\omega$ is a subexpression of $f$ and no subexpression of $f$ is more rapidly varying than $\omega$. Let $mrv(f)$ be the set of most rapidly varying subexpressions of $f$. We will write $mrv(f)\prec g$ (resp. $mrv(f)\asymp g$) if for all $\omega\in mrv(f)$, $\omega\prec g$ (resp. $\omega\asymp g$). Let $mrv(f_{1},\ldots,f_{n})$ be the set of $\omega$ such that $\omega$ is a subexpression of some $f_{i}$ and no subexpression of any $f_{j}$ is more rapidly varying than $\omega$, that is $mrv(f_{1},\ldots,f_{n})$ the $\max$ of $mrv(f_{1}),\ldots,mrv(f_{n})$ (as in Algorithm 3.12 of [2]).

To prove termination of our algorithm we use the notion of size of an exp-log expression defined in [2], section 3.4.1. For an exp-log expression $f$, let $S(f)$ be the set of subexpressions of $f$ defined by the following conditions.

1. (1)

If $f$ does not contain the variable $x$, then $S(f)=\emptyset$. 2. (2)

If $f=x$, then $S(f)=\{x\}$. 3. (3)

If $f=g+h$, $f=gh$, or $f=\frac{g}{h}$, then $S(f)=S(g)\cup S(h)$. 4. (4)

If $f=g^{c}$ then $S(f)=S(g)$. 5. (5)

If $f=\exp g$ or $f=\log g$ then $S(f)=\{f\}\cup S(g)$.

Then $Size(f)$ is defined as the cardinality of $S(f)$.

Let $\exp^{k}$ (resp. $\log^{k}$) denote $k$ times iterated exponential (resp. logarithm), and for $f\in E_{\infty}(x)$ let $f^{\uparrow k}$ (resp. $f^{\downarrow k}$) denote $f$ with $x$ replaced with $\exp^{k}(x)$ (resp. $\log^{k}(x)$). The following algorithm computes approximations of elements of a finite subset of $E_{\infty}(x)$ in terms of their most rapidly varying subexpression.

Algorithm 14**.**

*(MrvApprox)

Input: $a_{0},\ldots,a_{n}\in E_{\infty}(x)$ such that $\sum_{i=0}^{n}Size(a_{i})>0$.

Output: $\omega\in E_{\infty}(x)$, $b_{0},\ldots,b_{n}\in E_{\infty}(x)$, $e_{0},\ldots,e_{n}\in\mathbb{R}\cup\{\infty\}$, $d>0$, and $k\in\mathbb{Z}_{\geq 0}$ such that*

• •

$\omega>0$* and $\lim_{x\rightarrow\infty}\omega=0$,*

• •

$mrv(b_{0},\ldots,b_{n})\prec\omega$,

• •

if $a_{i}=0$ then $b_{i}=0$ and $e_{i}=\infty$,

• •

if $a_{i}\neq 0$ then $e_{i}\in\mathbb{R}$ and $\lim_{x\rightarrow\infty}\omega^{-(e_{i}+d)}(a_{i}^{\uparrow k}-b_{i}\omega^{e_{i}})=0$,

• •

$\sum_{i=0}^{n}Size(b_{i})<\sum_{i=0}^{n}Size(a_{i})$.

The algorithm proceeds in a very similar manner to the algorithm MrvLimit described in [2]. First, it finds the set $\Omega=mrv(a_{0},\ldots,a_{n})$. If $x\in\Omega$ the algorithm replaces $x$ with $\exp(x)$ in $a_{0},\ldots,a_{n}$ and recomputes $\Omega$ until $x\notin\Omega$. $k$ is the number of replacements performed in this step. Then the algorithm picks $\omega$ such that $\omega$ or $1/\omega$ belongs to $\Omega$, $\omega>0$, and $\lim_{x\rightarrow\infty}\omega=0$, and rewrites all elements of $\Omega$ in terms of $\omega$. If $a_{i}^{\uparrow k}$ contains a subexpression in the same comparability class as $\omega$, let $b_{i}\omega^{e_{i}}$ be the first term of $Series(a_{i}^{\uparrow k},\omega)$ (as in section 3.3.3 of [2]), and let $d_{i}>0$ be the difference between the exponents of $\omega$ in the second and in the first term of the series ($d_{i}=\infty$ if $a_{i}^{\uparrow k}=b_{i}\omega^{e_{i}}$). If $a_{i}^{\uparrow k}$ does not contain subexpressions in the same comparability class as $\omega$, then $b_{i}=a_{i}^{\uparrow k}$, $e_{i}=0$, and $d_{i}=\infty$. In both cases $mrv(b_{i})\prec\omega$. Pick $0<d<\min_{0\leq i\leq n}d_{i}$. Then $\omega^{-(e_{i}+d)}(a_{i}^{\uparrow k}-b_{i}\omega^{e_{i}})$ is either [math] or a power series in $\omega$ with positive exponents and coefficients in a lower comparability class than $\omega$, hence $\lim_{x\rightarrow\infty}\omega^{-(e_{i}+d)}(a_{i}^{\uparrow k}-b_{i}\omega^{e_{i}})=0$. Section 3.4.1 of [2] proves that $Size(b_{i})\leq Size(a_{i})$, with the strict inequality if $a_{i}^{\uparrow k}$ contains a subexpression in the same comparability class as $\omega$. This shows that the last requirement is satisfied.

3. Root continuity

To prove correctness of our algorithm we need a polynomial root continuity lemma that does not assume fixed degree of the polynomial. The lemma is very similar to Theorem 1 of [15], except that our version provides explicit bounds.

Lemma 15**.**

Let

[TABLE]

where $n\geq 1$ and $a_{n}\neq 0$. Let $\Gamma=\max_{1\leq i\leq n}\lvert r_{i}\rvert$ and $\Delta=\min_{r_{i}\neq r_{j}}\lvert r_{i}-r_{j}\rvert$ ($\Delta=\infty$ if $r_{1}=\ldots=r_{n}$).

Suppose that $m\geq n$, $0<\epsilon<\min(1,\frac{1}{\Gamma+1},\frac{\Delta}{2})$, and $0<\delta<\lvert a_{n}\rvert\frac{(1-\epsilon)\epsilon^{m+n}}{1-\epsilon^{m+1}}$.

Then for every

[TABLE]

such that $b_{m}\neq 0$, for $0\leq i\leq n$, $\lvert b_{i}-a_{i}\rvert<\delta$, and, for $n+1\leq i\leq m$, $\lvert b_{i}\rvert<\delta$, we have

[TABLE]

for $1\leq i\leq n$, $\lvert s_{i}-r_{i}\rvert<\epsilon$, and for $n+1\leq i\leq m$, $\lvert s_{i}\rvert>1/\epsilon$.

Proof.

Let $C_{0}=\{c\,:\,\lvert c\rvert=1/\epsilon\}$, $D_{0}=\{c\,:\,\lvert c\rvert\leq 1/\epsilon\}$, and, for $1\leq i\leq n$, let $C_{i}=\{c\,:\,\lvert c-r_{i}\rvert=\epsilon\}$, and $D_{i}=\{c\,:\,\lvert c-r_{i}\rvert\leq\epsilon\}$. Then, for $1\leq i,j\leq n$, $D_{i}$ and $D_{j}$ are either identical or disjoint, $D_{i}$ is contained in the interior of $D_{0}$, $D_{i}$ contains exactly one of the distinct roots of $p$, and $D_{0}$ contains all roots of $p$. If $c\in C_{i}$ for some $0\leq i\leq n$, then $\lvert p(c)\rvert\geq\lvert a_{n}\rvert\epsilon^{n}$ and

[TABLE]

We have

[TABLE]

Hence $\lvert q(c)-p(c)\rvert<\lvert p(c)\rvert$. By Rouche’s theorem, for $0\leq i\leq n$, the number of roots of $q$ in $D_{i}$ equals the number of roots of $p$ in $D_{i}$, which concludes the proof. ∎

4. The main algorithm

Let $E_{\infty}^{\mathbb{C}}(x)=\{u+\imath v\,:\,u,v\in E_{\infty}(x)\}$. This section presents the main algorithm computing asymptotic approximations of roots of polynomials $P\in E_{\infty}^{\mathbb{C}}(x)[y]$.

Let us first describe a straightforward generalization of Algorithm 14 to inputs in $E_{\infty}^{\mathbb{C}}(x)$. We extend $mrv$ and $Size$ to $E_{\infty}^{\mathbb{C}}(x)$ by defining $mrv(u_{1}+\imath v_{1},\ldots,u_{n}+\imath v_{n})=mrv(u_{1},v_{1},\ldots,u_{n},v_{n})$ and $Size(u+\imath v)=Size(u)+Size(v)$. We say that $u+\imath v\prec\omega$ (resp. $u+\imath v\asymp\omega$) if $\lvert u+\imath v\rvert=(u^{2}+v^{2})^{1/2}\prec\omega$ (resp. $\lvert u+\imath v\rvert\asymp\omega$). If $a=u+\imath v$ then $a^{\uparrow k}:=u^{\uparrow k}+\imath v^{\uparrow k}$ and $a^{\downarrow k}:=u^{\downarrow k}+\imath v^{\downarrow k}$.

Algorithm 16**.**

*(MrvApproxC)

Input: $a_{0},\ldots,a_{n}\in E_{\infty}^{\mathbb{C}}(x)$ such that $\sum_{i=0}^{n}Size(a_{i})>0$.

Output: $\omega\in E_{\infty}(x)$, $b_{0},\ldots,b_{n}\in E_{\infty}^{\mathbb{C}}(x)$, $e_{0},\ldots,e_{n}\in\mathbb{R}\cup\{\infty\}$, $d>0$, and $k\in\mathbb{Z}_{\geq 0}$ such that*

• •

$\omega>0$* and $\lim_{x\rightarrow\infty}\omega=0$,*

• •

$mrv(b_{0},\ldots,b_{n})\prec\omega$,

• •

if $a_{i}=0$ then $b_{i}=0$ and $e_{i}=\infty$,

• •

if $a_{i}\neq 0$ then $e_{i}\in\mathbb{R}$ and $\lim_{x\rightarrow\infty}\omega^{-(e_{i}+d)}(a_{i}^{\uparrow k}-b_{i}\omega^{e_{i}})=0$,

• •

$\sum_{i=0}^{n}Size(b_{i})<\sum_{i=0}^{n}Size(a_{i})$.

1. (1)

Let $a_{i}=u_{i}+\imath v_{i}$. Call Algorithm 14 with $u_{0},v_{0},\ldots,u_{n},v_{n}$ as input, obtaining $\omega\in E_{\infty}(x)$, $\hat{u}_{0},\hat{v}_{0},\ldots,\hat{u}_{n},\hat{v}_{n}\in E_{\infty}(x)$, $e_{u,0},e_{v,0},\ldots,e_{u,n}e_{v,n}\in\mathbb{R}$, $\hat{d}>0$, and $k\in\mathbb{Z}_{\geq 0}$. 2. (2)

For $0\leq i\leq n$, put $e_{i}:=\min(e_{u,i},e_{v,i})$ and

[TABLE] 3. (3)

Pick $d$ such that $0<d<\hat{d}$ and $d<\lvert e_{u,i}-e_{v,i}\rvert$ for all $i$ such that $e_{u,i}\neq e_{v,i}$. 4. (4)

Return $\omega$, $b_{0},\ldots,b_{n}$, $e_{0},\ldots,e_{n}$, $d$, and $k$.

Proof.

To prove that the output of Algorithm 16 satisfies the required conditions we need to prove that if $a_{i}\neq 0$ then

[TABLE]

The other conditions follow directly from the definitions and the properties of the output of Algorithm 14.

Suppose that $e_{u,i}<e_{v,i}$. Then $e_{i}=e_{u,i}$ and

[TABLE]

We have

[TABLE]

and

[TABLE]

since $(e_{v,i}-e_{u,i})+(\hat{d}-d)>0$, $(e_{v,i}-e_{u,i})-d>0$, and $\hat{v}_{i}\prec\omega$. Cases $e_{u,i}=e_{v,i}$ and $e_{u,i}>e_{v,i}$ can be proven in a similar manner. ∎

Let $P(x,y)=a_{n}(x)y^{n}+\ldots+a_{0}(x)\in E_{\infty}^{\mathbb{C}}(x)[y]$. W.l.o.g. we may assume that $a_{n}$ and $a_{0}$ are not identically zero.

Suppose that $\sum_{i=0}^{n}Size(a_{i})>0$ i.e. $P(x,y)$ depends on $x$. Let $\omega\in E_{\infty}(x)$, $b_{0},\ldots,b_{n}\in E_{\infty}^{\mathbb{C}}(x)$, $e_{0},\ldots,e_{n}\in\mathbb{R}\cup\{\infty\}$, $d>0$, and $k\in\mathbb{Z}_{\geq 0}$ be the output of Algorithm 16 for $a_{0},\ldots,a_{n}$, and let $Q(x,y)=a_{n}^{\uparrow k}(x)y^{n}+\ldots+a_{0}^{\uparrow k}(x)$.

Let $K$ be a Hardy field containing germs at infinity of exp-log functions defined near infinity, such that $K[\imath]$ is algebraically closed. Let $\alpha\in K[\imath]$ be a root of $Q$. Since $\lvert\alpha\rvert\in K$, the limit $\lim_{x\rightarrow\infty}\frac{\log\lvert\alpha\rvert}{\log\lvert\omega\rvert}=\gamma$ exists.

Claim 17*.*

$\gamma\in\mathbb{R}$.

Proof.

Suppose that $\lvert\gamma\rvert=\infty$. Then

[TABLE]

Since $\omega\prec\lvert\alpha\rvert$, $\frac{a_{n-i}^{\uparrow k}}{a_{n}^{\uparrow k}}\prec\lvert\alpha\rvert$, and so either $\lim_{x\rightarrow\infty}\lvert\alpha\rvert=\infty$ and

[TABLE]

or $\lim_{x\rightarrow\infty}\lvert\alpha\rvert=0$ and

[TABLE]

Both cases contradict equation (4.1). ∎

Let $I=\{i\,:\,0\leq i\leq n\wedge a_{i}\neq 0\}$. If $i\in I$, put

[TABLE]

Then $a_{i}^{\uparrow k}=b_{i}\omega^{e_{i}}+c_{i}\omega^{e_{i}+d}$ and $\lim_{x\rightarrow\infty}c_{i}=0$. Put $\beta:=\frac{\alpha}{\omega^{\gamma}}$. Then

[TABLE]

and hence $\beta\prec\omega$. We have

[TABLE]

Let $\mu=\min_{i\in I}e_{i}+\gamma i$, let $J=\{i\in I\,:\,e_{i}+\gamma i=\mu\}$, and let

[TABLE]

Then

[TABLE]

As $x$ tends to infinity, all terms in the last two sums tend to zero, hence

[TABLE]

Since $\beta\prec\omega$ and $\mu-\nu<0$, the cardinality of $J$ must be at least $2$. Consider the subset $A=\{(i,e_{i})\,:\,i\in I\}$ of $\mathbb{R}^{2}$ with coordinates denoted $(x_{1},x_{2})$. Then the line $x_{2}=-\gamma x_{1}+\mu$ passes through the points $\{(i,e_{i})\,:\,i\in J\}$ and all the other points of $A$ lie above this line. This means that $-\gamma$ is the slope of one of the segments that form the lower part of the boundary of the convex hull of $A$.

Let $R_{\gamma}(x,z)=\sum_{i\in J}b_{i}z^{i}$ and $Q_{\gamma}(x,z)=\omega^{-\mu}Q(x,z\omega^{\gamma})$. We have

[TABLE]

Let $\rho_{1},\ldots,\rho_{t}\in K[\imath]\setminus\{0\}$ be the nonzero roots of $R_{\gamma}(x,z)$ listed with multiplicities.

Claim 18*.*

There exist roots $\beta_{1},\ldots,\beta_{t}\in K[\imath]$ of $Q_{\gamma}(x,z)$ such that, for sufficiently large $x$, for $1\leq j\leq t$, $\lvert\beta_{j}-\rho_{j}\rvert<\omega^{\eta}$, where $\eta=\frac{\nu-\mu}{4(n+1)}>0$.

Proof.

Let $\Gamma(x)$ be the maximum of absolute values of roots of $R_{\gamma}(x,z)$ and let $\Delta(x)$ be the minimum distance between two distinct roots of $R_{\gamma}(x,z)$ ($\Delta(x)=\infty$ if all roots of $R_{\gamma}(x,z)$ are equal). Put $\delta(x)=\omega(x)^{(\nu-\mu)/2}$ and $\epsilon(x)=\delta(x)^{1/(2n+2)}$. For sufficiently large $x$, $R_{\gamma}(x,z)$ has a fixed number of distinct roots in $z$, equal to its number of distinct roots in $K[\imath]$. $\Gamma(x)$ can be bounded from above by a rational function in absolute values of $b_{i}$, for $i\in J$, and, since the coefficients in $z$ of $R_{1}(x,z)/g.c.d.(R_{1}(x,z),\frac{\partial}{\partial z}R_{1}(x,z))$ are rational functions of $b_{i}$, $\Delta(x)$ can be bounded from below by an expression constructed from $b_{i}$ using rational operations, square roots and absolute value (see e.g. [5], Theorem 5). Since $mrv(b_{i})\prec\omega$, for sufficiently large $x$, $0<\epsilon<\min(1,\frac{1}{\Gamma+1},\frac{\Delta}{2})$. Let $s$ be the degree of $R_{\gamma}$ in $z$. For sufficiently large $x$,

[TABLE]

The coefficients at $z^{i}$, for $0\leq i\leq n$, of $Q_{\gamma}(x,z)-R_{\gamma}(x,z)$ have the form $\omega^{\nu-\mu}\xi_{i}$ and $\lim_{x\rightarrow\infty}\xi_{i}=0$, hence, for sufficiently large $x$, the absolute value of each of these coefficients is less than $\delta$. By Lemma 15, there exist roots $\beta_{1},\ldots,\beta_{t}\in K[\imath]$ of $Q_{\gamma}(x,z)$ such that, for sufficiently large $x$, for $1\leq j\leq t$, $\lvert\beta_{j}-\rho_{j}\rvert<\epsilon=\omega^{\eta}$. ∎

Fix $1\leq j\leq t$ and let $\alpha_{j}=\beta_{j}\omega^{\gamma}$. Then $Q(x,\alpha_{j})=\omega^{\mu}Q_{\gamma}(x,\beta_{j})=0$. Suppose that $f_{1},\ldots,f_{m}$ form an $m$-term asymptotic approximation* of *$\rho_{j}$ such that, for $1\leq i\leq m$, $mrv(f_{i})\prec\omega$. For $1\leq i\leq m$, put $g_{i}=f_{i}\omega^{\gamma}$. We have

[TABLE]

For sufficiently large $x$, $\lvert\beta_{j}-\rho_{j}\rvert<\omega^{\eta}$. Since $f_{m}\prec\omega$, $\lim_{x\rightarrow\infty}\frac{\omega^{\eta}}{f_{m}}=0$. Hence, $\lim_{x\rightarrow\infty}\frac{\alpha_{j}-\sum_{i=1}^{m}g_{i}}{g_{m}}=0$, and so $g_{1},\ldots,g_{m}$ form an $m$-term asymptotic approximation* of *$\alpha_{j}$. Since for any $f\in K[\imath]$

[TABLE]

$g_{1}^{\downarrow k},\ldots,g_{m}^{\downarrow k}$ form an $m$-term asymptotic approximation of the root $\alpha_{j}(\log^{k}(x))$ of $P$.

Let us now consider the case where we have found an exact solution $\rho_{j}\in E_{\infty}^{\mathbb{C}}(x)$ of $R_{\gamma}(x,z)$. To simplify the description of the case let us make the following rather technical definition.

Definition 19**.**

Let $\omega\in E_{\infty}(x)$, $b_{0},\ldots,b_{n}\in E_{\infty}^{\mathbb{C}}(x)$, $e_{0},\ldots,e_{n}\in\mathbb{R}\cup\{\infty\}$, $d>0$, and $k\in\mathbb{Z}_{\geq 0}$ be the result of applying Algorithm 16 to $a_{0},\ldots,a_{n}$. We will call a root $\alpha(\log^{k}(x))$ of $P$ asymptotically small if $\lim_{x\rightarrow\infty}\alpha=0$ and $\alpha\asymp\omega$ (in other words, $\alpha=\beta\omega^{\gamma}$ with $\gamma>0$).

Suppose that we have found an exact solution $\rho_{j}\in E_{\infty}^{\mathbb{C}}(x)$ of $R_{\gamma}(x,z)$ of multiplicity $\sigma_{j}$. Then there exist exactly $\sigma_{j}$ roots $\beta_{1,j},\ldots,\beta_{\sigma_{j},j}\in K[\imath]$ of $Q_{\gamma}(x,z)$ such that, for sufficiently large $x$, for $1\leq\iota\leq\sigma_{j}$, $\lvert\beta_{\iota,j}-\rho_{j}\rvert<\omega^{\eta}$. Hence, there exist exactly $\sigma_{j}$ roots $\alpha_{\iota,j}=\beta_{\iota,j}\omega^{\gamma}$ of $Q(x,y)$ such that, for sufficiently large $x$, for $1\leq\iota\leq\sigma_{j}$, $\lvert\alpha_{\iota,j}-\rho_{j}\omega^{\gamma}\rvert<\omega^{\gamma+\eta}$. The mapping $\varphi:\zeta\rightarrow\rho_{j}\omega^{\gamma}+\omega^{\gamma}\zeta$ is a bijection between the roots $\zeta$ of $Q(x,\rho_{j}\omega^{\gamma}+\omega^{\gamma}y)$ such that, for sufficiently large $x$, $\lvert\zeta\rvert<\omega^{\eta}$ and the roots $\varphi(\zeta)$ of $Q(x,y)$ such that $\lvert\varphi(\zeta)-\rho_{j}\omega^{\gamma}\rvert<\omega^{\gamma+\eta}$. Since $\nu>\mu$ can be chosen arbitrarily close to $\mu$, $\eta$ can be arbitrarily small. Therefore the mapping $\xi\rightarrow(\rho_{j}\omega^{\gamma})^{\downarrow k}+(\omega^{\gamma})^{\downarrow k}\xi$ is a bijection between the roots $\xi$ of $P(x,(\rho_{j}\omega^{\gamma})^{\downarrow k}+(\omega^{\gamma})^{\downarrow k}y)$ that are identically zero or asymptotically small and the roots $\alpha_{1,j}(\log^{k}(x)),\ldots,\alpha_{\sigma_{j},j}(\log^{k}(x))$ of $P(x,y)$.

The above discussion suggests the following procedure for finding asymptotic approximations of roots of $P$. Use Algorithm 16 for $a_{0},\ldots,a_{n}$, to find $\omega\in E_{\infty}(x)$, $b_{0},\ldots,b_{n}\in E_{\infty}^{\mathbb{C}}(x)$, $e_{0},\ldots,e_{n}\in\mathbb{R}\cup\{\infty\}$, $d>0$, and $k\in\mathbb{Z}_{\geq 0}$. Compute the values of $\gamma$ such that $-\gamma$ is the slope of one of the segments that form the lower part of the boundary of the convex hull of $A=\{(i,e_{i})\,:\,i\in I\}$. For each $\gamma$ find $R_{\gamma}(x,z)$ and call the procedure recursively to find asymptotic approximations of roots of $R_{\gamma}$. Finally, obtain asymptotic approximations of roots of $P$ by multiplying the terms of asymptotic approximations of roots of $R_{\gamma}$ by $\omega^{\gamma}$ and replacing $x$ with $\log^{k}(x)$. The following algorithm formalizes this procedure, handles the base case, and the case where we get an exact solution with less than the requested $m$ terms.

Notation 20*.*

We use the notation $\sqcup$ for joining lists, that is

[TABLE]

For a list $F=(f_{1},\ldots,f_{l})$ of expressions in $E_{\infty}^{\mathbb{C}}(x)$ let

[TABLE]

let

[TABLE]

and, for $g\in E_{\infty}^{\mathbb{C}}(x)$, let

[TABLE]

Algorithm 21**.**

*(AsymptoticSolutions)

Input: $P(x,y)=a_{n}(x)y^{n}+\ldots+a_{0}(x)\in E_{\infty}^{\mathbb{C}}(x)[y]$ with $a_{n}\neq 0$ and $a_{0}\neq 0$, $m\in\mathbb{Z}_{>0},$ $\operatorname{sflag}\in\{\operatorname{true},\operatorname{false}\}$

Output: $((F_{1},\sigma_{1}),\ldots,(F_{t},\sigma_{t}))$ such that*

• •

for $1\leq i\leq t$,* $F_{i}=(f_{i,1},\ldots,f_{i,m_{i}})$ is an $m_{i}$-term asymptotic approximation of $\sigma_{i}$ roots of $P(x,y)$ in $y$ (counted with multiplicities),*

• •

$f_{i,1},\ldots,f_{i,m_{i}}\in E_{\infty}^{\mathbb{C}}(x)$,

• •

either $m_{i}=m$ or $m_{i}<m$ and $f_{i,1}+\ldots+f_{i,m_{i}}$ is an exact root of $P$ of multiplicity $\sigma_{i}$,* ***

• •

if $\operatorname{sflag}=\operatorname{false}$ then $\sigma_{1}+\ldots+\sigma_{t}=n$ and $F_{1},\ldots,F_{t}$ are asymptotic approximations of all complex roots of $P(x,y)$,

• •

if $\operatorname{sflag}=\operatorname{true}$ then $F_{1},\ldots,F_{t}$ are asymptotic approximations of all asymptotically small complex roots of $P(x,y)$.

1. (1)

If $P$ does not depend on $x$ then

1. (a)

if $\operatorname{sflag}=\operatorname{true}$ return $()$, 2. (b)

let $r_{1},\ldots,r_{t}\in\mathbb{C}$ be the distinct roots of $P$, 3. (c)

for $1\leq i\leq t$,* let $\sigma_{i}$ be the multiplicity of $r_{i}$,* 4. (d)

return $(((r_{1}),\sigma_{1}),\ldots,((r_{t}),\sigma_{t}))$. 2. (2)

Apply Algorithm 16 to $a_{0},\ldots,a_{n}$, obtaining $\omega\in E_{\infty}(x)$,

[TABLE]

$e_{0},\ldots,e_{n}\in\mathbb{R}\cup\{\infty\}$, $d>0$, and $k\in\mathbb{Z}_{\geq 0}$. 3. (3)

Let $I=\{i\,:\,0\leq i\leq n\wedge b_{i}\neq 0\}$* and $A=\{(i,e_{i})\>:\>i\in I\}$. Compute $\gamma_{1},\ldots,\gamma_{l}$ such that the lower part of the boundary of the convex hull of $A$ consists of segments with slopes $-\gamma_{1},\ldots,-\gamma_{l}$.* 4. (4)

Set $\mathcal{R}=()$. 5. (5)

For $1\leq j\leq l$ do:

1. (a)

if $\operatorname{sflag}=\operatorname{true}$ and $\gamma_{j}\leq 0$, continue the loop with the next $j$, 2. (b)

compute $\mu=\min_{i\in I}e_{i}+\gamma_{j}i$ , $J=\{i\in I\,:\,e_{i}+\gamma_{j}i=\mu\}$, 3. (c)

let $\lambda=\min J$ and let $R(x,z)=\sum_{i\in J}b_{i}z^{i-\lambda}$, 4. (d)

compute

[TABLE]

where $F_{\iota}=(f_{\iota,1},\ldots,f_{\iota,m_{\iota}})$, for $1\leq\iota\leq t$, 5. (e)

For $1\leq\iota\leq t$ do:

1. (i)

put $G=(\omega^{\gamma_{j}}F)^{\downarrow k}$, 2. (ii)

if $m_{\iota}=m$, set $\mathcal{R}=\mathcal{R}\sqcup((G,\sigma_{\iota}))$ and continue the loop with the next $\iota$, 3. (iii)

if $m_{\iota}<m$, put $r=\sum G$, 4. (iv)

compute $P(x,r+(\omega^{\gamma_{j}})^{\downarrow k}y)=a_{r,n}(x)y^{n}+\ldots+a_{r,0}(x)$, with $a_{r,i}\in E_{\infty}^{\mathbb{C}}(x)$, and let $\lambda=\min\{i\,:\,a_{r,i}\neq 0\}$, 5. (v)

if $\lambda>0$, set $\mathcal{R}=\mathcal{R}\sqcup((G,\lambda))$, and if $\lambda=\sigma_{\iota}$, continue the loop with the next $\iota$, 6. (vi)

put $P_{r}(x,y)=a_{r,n}(x)y^{n-\lambda}+\ldots+a_{r,\lambda}(x)$, 7. (vii)

compute

[TABLE] 8. (viii)

for $1\leq\kappa\leq s$, 9. (ix)

set

[TABLE] 6. (6)

Return $\mathcal{R}$.

Proof.

For a proof of termination of the algorithm, put $\Sigma=\sum_{i=0}^{n}Size(a_{i})$ and let us use pairs $(m,\Sigma)\in\mathbb{Z}_{>0}\times\mathbb{Z}_{\geq 0}$ as a metric. Note that $\mathbb{Z}_{>0}\times\mathbb{Z}_{\geq 0}$ with the lexicographic order does not admit infinite strictly decreasing sequences. The recursive calls to Algorithm 21 in step $5(d)$ have the same value of $m$ and a strictly lower value of $\Sigma$, because according to the specification of Algorithm 16, $\sum_{i=0}^{n}Size(b_{i})<\sum_{i=0}^{n}Size(a_{i})$. The recursive calls to Algorithm 21 in step $5(e)(vii)$ have a strictly lower value of $m$. This shows that the algorithm terminates.

Correctness of the algorithm follows from the discussion earlier in this section. ∎

Example 22**.**

Find one-term asymptotic approximations of roots of

[TABLE]

Applying Algorithm 16 to the coefficients of $P$ we get $\omega=\exp(-x)$,

[TABLE]

$d=\infty$, and $k=0$ (the coefficients and exponents can be read off $P$ written in terms of $\omega$: $y^{5}-\omega^{-1}y^{4}+x\omega^{-\pi}y^{3}+\log(x)y-x^{2}$; in this case $a_{i}^{\uparrow k}=b_{i}\omega^{e_{i}}$ hence $d=\infty$). The set $A$ and the lower part of the boundary of the convex hull of $A$ are shown in Figure 4.1. We obtain $\gamma_{1}=\frac{\pi}{3}$ and $\gamma_{2}=-\frac{\pi}{2}$.

For $\gamma_{1}$ we get $R_{1}=xy^{3}-x^{2}$. In the recursive call to Algorithm 21 applying Algorithm 16 to the coefficients of $R_{1}$ yields $\omega_{1}=\exp(-x)$,

[TABLE]

$d_{1}=\infty$, and $k_{1}=1$. The lower part of the boundary of the convex hull of $A_{1}$ consists of one segment and we get $\gamma_{1,1}=-\frac{1}{3}$ and the corresponding polynomial $R_{1,1}=y^{3}-1$. The recursive call to Algorithm 21 returns simple roots $1,\frac{-1-\imath\sqrt{3}}{2},\frac{-1+\imath\sqrt{3}}{2}$ of $R_{1,1}$. We have $(\omega_{1}^{\gamma_{1,1}})^{\downarrow 1}=\sqrt[3]{x}$ hence Algorithm 21 for $R_{1}$ returns

[TABLE]

all with multiplicity $1$. Since $(\omega^{\gamma_{1}})^{\downarrow 0}=\exp(-\frac{\pi}{3}x)$, we add

[TABLE]

to $\mathcal{R}$.

For $\gamma_{2}$ we get $R_{2}=y^{2}+x$. In the recursive call to Algorithm 21 applying Algorithm 16 to the coefficients of $R_{2}$ yields $\omega_{2}=\exp(-x)$, $(b_{2,0},b_{2,1},b_{2,2})=(1,0,1)$, $(e_{2,0},e_{2,1},e_{2,2})=(-1,\infty,0)$, $d_{2}=\infty$, and $k_{2}=1$. The lower part of the boundary of the convex hull of $A_{2}$ consists of one segment and we get $\gamma_{2,1}=-\frac{1}{2}$ and the corresponding polynomial $R_{2,1}=y^{2}+1$. The recursive call to Algorithm 21 returns simple roots $-\imath,\imath$ of $R_{2,1}$. We have $(\omega_{2}^{\gamma_{2,1}})^{\downarrow 1}=\sqrt{x}$ hence Algorithm 21 for $R_{2}$ returns $-\imath\sqrt{x},\imath\sqrt{x}$, both with multiplicity $1$. Since $(\omega^{\gamma_{2}})^{\downarrow 0}=\exp(\frac{\pi}{2}x)$, we add $-\imath\sqrt{x}\exp(\frac{\pi}{2}x),\imath\sqrt{x}\exp(\frac{\pi}{2}x)$ to $\mathcal{R}$.

Finally, the algorithm returns one-term asymptotic approximations

[TABLE]

all with multiplicity $1$.

5. Real roots

Imaginary part of an non-real root may be asymptotically smaller than the real part, hence real asymptotic approximations can correspond to non-real roots.

Example 23**.**

Let $P(x,y)=(y^{2}-x\exp(x)y+\exp(2x))^{2}+1$. Three-term asymptotic approximations of roots of $P$ in $y$ computed with Algorithm 21 are $(x^{-1}+x^{-3}+2x^{-5})\exp(x)$ and $(x-x^{-1}-x^{-3})\exp(x)$, both with multiplicity two. The approximations are real-valued, yet $P$ clearly does not have real roots. Computing more terms only adds real-valued terms of the form $ax^{-n}$ in the coefficient of $\exp(x)$. Imaginary parts would show up only in a transfinite series representation, since the imaginary parts are asymptotically smaller than $x^{-n}\exp(x)$ for any $n$. For this low degree polynomial $P$ we can compute asymptotic approximations of imaginary parts of roots of $P$, by computing $R=res_{z}(P(x,y+z),P(x,z))$. The roots of this polynomial of degree $16$ in $y$ with $32$ terms, are differences of pairs of roots of $P$. In particular, four of the roots are equal to the imaginary parts of roots of $P$, multiplied by two. One-term asymptotic approximations of roots of $R$ computed with Algorithm 21 include two purely imaginary-valued expressions, $-2\iota x^{-1}\exp(-x)$ and $2\iota x^{-1}\exp(-x)$, both of multiplicity two. This shows that, indeed, the imaginary parts of roots of $P$ are asymptotically smaller than $x^{-n}\exp(x)$ for any $n$.

In this section we provide a method for deciding which asymptotic approximations correspond to real roots.

Let $P(x,y)=a_{n}(x)y^{n}+\ldots+a_{0}(x)\in E_{\infty}(x)[y]$ with $a_{n}\neq 0$ and $a_{0}\neq 0$, and let $((F_{1},\sigma_{1}),\ldots,(F_{t},\sigma_{t}))$ be the output of Algorithm 21 for $P$, with some $m>0$ and $sflag=false$. Note that here we assume that the coefficients of $P$ are real-valued. For $1\leq i\leq t$, $F_{i}=(f_{i,1},\ldots,f_{i,m_{i}})$. First, let us note the easy cases. If any of $f_{i,j}$ are not real-valued, then $F_{i}$ does not correspond to a real root. If all $f_{i,j}$ are real-valued and $\sigma_{i}=1$ then $F_{i}$ corresponds to a real root. If $m_{i}<m$ then $f_{i,1}+\ldots+f_{i,m_{i}}$ is equal to the exact root, hence it is evident whether the root is real valued.

The hard case is when there are real-valued asymptotic approximations with $\sigma_{i}>1$ and $m_{i}=m$. We will find the number of distinct real roots corresponding to each real-valued asymptotic approximation with multiplicity higher than one. Assume that, possibly after reordering, the real-valued asymptotic approximations are $F_{i}$, for $1\leq i\leq s$. Let $r_{i}=f_{i,1}+\ldots+f_{i,m_{i}}$, for $1\leq i\leq s$. The algorithm MrvLimit of [2] contains a subprocedure which computes the sign of exp-log expressions near infinity. Hence, we can reorder the approximations so that, for sufficiently large $x$, $r_{1}<\ldots<r_{s}$. Put $h_{0}=-\infty$, $h_{i}=\frac{r_{i}+r_{i+1}}{2}$, for $1\leq i<s$, and $h_{s}=\infty$.

Lemma 24**.**

If $F_{i}$ corresponds to a real root $\alpha$ of $P$ then, for sufficiently large $x$, $h_{i-1}<\alpha<h_{i}$.

To prove the lemma we will use the following claim.

Claim 25*.*

If $F=(f_{1},\ldots,f_{m_{f}})$ and $G=(g_{1},\ldots,g_{m_{g}})$ are asymptotic approximations returned by Algorithm 21 for $P$ and there is $l\leq\min(m_{f},m_{g})$ such that, for all $1\leq j<l$, $f_{j}=g_{j}$ and $f_{l}\neq g_{l}$, then $\lim_{x\rightarrow\infty}\frac{f_{l}}{g_{l}}\neq 1$.

Proof.

The claim is true when $P$ does not depend on $x$, hence, by induction, we may assume that the claim is true for the recursive calls to Algorithm 21. If $F$ and $G$ were computed in different iterations of the loop in step $5$ then $l=1$, $f_{1}=(\omega^{\gamma_{f}}c_{f})^{\downarrow k}$, $g_{1}=(\omega^{\gamma_{g}}c_{g})^{\downarrow k}$, $mrv(c_{f},c_{g})\prec\omega$, and $\gamma_{f}\neq\gamma_{g}$, hence $\lim_{x\rightarrow\infty}\frac{\lvert f_{l}\rvert}{\lvert g_{l}\rvert}$ is either [math] or $\infty$. If $F$ and $G$ were computed in the same iteration of the loop in step $5$, but in different iterations of the loop in step $5(e)$ then the claim is true by the inductive hypothesis applied to $R(x,z)$. Finally, if $F$ and $G$ were computed in the same iteration of the loop in step $5(e)$ then $l>m_{\iota}$, and hence neither $F$ nor $G$ was added in step $5(e)(v)$. Therefore, the claim is true by the inductive hypothesis applied to $P_{r}(x,y)$. ∎

Let us now prove Lemma 24.

Proof.

If $m_{i}<m$ then $\alpha=r_{i}$ and $h_{i-1}<r_{i}<h_{i}$, because, for sufficiently large $x$, $r_{1}<\ldots<r_{s}$. Hence we can assume that $m_{i}=m$. Let us prove that $h_{i-1}<\alpha$. If $i=1$, then $h_{i-1}=-\infty$ and the inequality is true. Let $l$ be such that for all $1\leq j<l$ $f_{i-1,j}=f_{i,j}$ and $f_{i-1,l}\neq f_{i,l}$. We define $f_{i-1,m_{i-1}+1}=0$, so that such $l$ always exist. Note that, for sufficiently large $x$, $f_{i-1,l}<f_{i,l}$, because $r_{i-1}<r_{i}$. If $\lim_{x\rightarrow\infty}\frac{\lvert f_{i-1,l}\rvert}{\lvert f_{i,l}\rvert}\leq 1$ put $g=f_{i,l}$ else put $g=f_{i-1,l}$. We have

[TABLE]

Since $F_{i}$ is an asymptotic approximation of $\alpha$, we have

[TABLE]

and, for $j>l$, we have $\lim_{x\rightarrow\infty}\frac{f_{i,j}}{g}=0$ and $\lim_{x\rightarrow\infty}\frac{f_{i-1,j}}{g}=0$, therefore

[TABLE]

If $f_{i-1,l}=0$, then $g=f_{i,l}$, for sufficiently large $x$, $f_{i,l}>0$, and $\lim_{x\rightarrow\infty}\frac{f_{i,l}-f_{i-1,l}}{2g}=\frac{1}{2}$, hence for sufficiently large $x$, $\alpha-h_{i-1}>0$.

If $f_{i-1,l}\neq 0$, then the assumptions of Claim 25 are satisfied, and hence $\lim_{x\rightarrow\infty}\frac{f_{i-1,l}}{f_{i,l}}\neq 1$.

Suppose that $g=f_{i,l}$. Then

[TABLE]

Since for sufficiently large $x$, $f_{i,l}-f_{i-1,l}>0$, $\lim_{x\rightarrow\infty}\frac{\lvert f_{i-1,l}\rvert}{\lvert f_{i,l}\rvert}\leq 1$, and $\lim_{x\rightarrow\infty}\frac{f_{i-1,l}}{f_{i,l}}\neq 1$, hence, for sufficiently large $x$, $f_{i,l}>0$. Therefore,

[TABLE]

which shows that, for sufficiently large $x$, $h_{i-1}<\alpha$.

Now suppose that $g=f_{i-1,l}$. Then

[TABLE]

Since for sufficiently large $x$, $f_{i,l}-f_{i-1,l}>0$, and $\lim_{x\rightarrow\infty}\frac{\lvert f_{i-1,l}\rvert}{\lvert f_{i,l}\rvert}>1$, hence, for sufficiently large $x$, $f_{i-1,l}<0$. Therefore,

[TABLE]

which shows that, for sufficiently large $x$, $h_{i-1}<\alpha$. The proof that, for sufficiently large $x$, $\alpha<h_{i}$ is similar. ∎

Let $P_{1},\ldots,P_{k}$ be the Sturm sequence of $P$ in $y$ over $E_{\infty}(x)$ (that is coefficients that are identically zero near infinity are set to zero). For $1\leq i\leq s-1$ and $1\leq j\leq k$, $P_{j}(x,h_{i})\in E_{\infty}(x)$, hence it has a constant sign $\theta_{i,j}$ near infinity, and we can compute $\theta_{i,j}$ using a subprocedure of MrvLimit. Let $c_{j}y^{n_{j}}$ be the leading term of $P_{j}$, for $1\leq j\leq k$, let $\theta_{0,j}$ be the sign near infinity of $(-1)^{n_{j}}c_{j}$, and let $\theta_{s,j}$ be the sign near infinity of $c_{j}$. For $1\leq i\leq s-1$, let $\nu_{i}$ be the number of sign changes in the sequence $\Theta_{i}=(\theta_{i,1},\ldots,\theta_{i,k})$.

Criterion 26**.**

The number of distinct real roots of $P$ corresponding to the asymptotic approximation $F_{i}$ is equal to $\nu_{i-1}-\nu_{i}$.

Correctness of the criterion follows from Lemma 24 and Sturm’s theorem.

Example 27**.**

As in Example 23, let $P(x,y)=(y^{2}-x\exp(x)y+\exp(2x))^{2}+1$. One-term asymptotic approximations of roots of $P$ in $y$ computed with Algorithm 21 are $r_{1}=x^{-1}\exp(x)$ and $r_{2}=x\exp(x)$, both with multiplicity two. The Sturm sequence of $P$ in $y$ is

[TABLE]

We have $\Theta_{0}=(1,-1,1,1,1)$ and $\Theta_{2}=(1,1,1,-1,1)$. Since $\nu_{0}=\nu_{2}=2$, $P$ has no real roots near infinity (and we do not need to compute $\nu_{1}$, since it must equal $2$ as well).

Let $Q(x,y)=(y^{2}-x\exp(x)y+\exp(2x))^{2}-1$. One-term asymptotic approximations of roots of $Q$ in $y$ computed with Algorithm 21 are the same as for $P$. The Sturm sequence of $Q$ in $y$ is

[TABLE]

We have $\Theta_{0}=(1,-1,1,-1,1)$ and $\Theta_{2}=(1,1,1,1,1)$. Since $\nu_{0}=4$ and $\nu_{2}=0$, $P$ has four distinct real roots near infinity. This again is sufficient to tell that $r_{1}$ or $r_{2}$ correspond to two real root each. And indeed, if we substitute $h_{1}=(r_{1}+r_{2})/2=(x^{-1}+x)\exp(x)$ into the Sturm sequence and compute the signs near infinity we get $\Theta_{1}=(1,-1,-1,1,1)$ and $\nu_{1}=2$.

6. Implementation and experimental results

We have implemented AsymptoticSolutions as a part of the Mathematica system. The implementation has been done in Wolfram Language, using elements of the MrvLimit algorithm, which is implemented partly in the C source code of Mathematica and partly in Wolfram Language. The experiments have been run on a laptop computer with a $2.7$ GHz Intel Core i7-4800MQ processor and $10$ GB of RAM assigned to the Linux virtual machine.

Example 28**.**

We use eight exp-log expressions from examples in [2] as polynomial coefficients.

[TABLE]

Let $P_{n}(x,y)=\sum_{i=0}^{n}a_{i}y^{i}$. We have run the examples with $n$ ranging from $2$ to $7$ and with varying number $m$ of requested terms. The results are given in Table 1. For each value of $m$ the row Time gives the computation time in seconds, the row Iter gives the number of calls to Algorithm 21, and the row LC gives the total leaf count of the returned expressions.

We can observe that for a fixed polynomial the number of recursive calls is close to linear in the number of additional terms requested. Increasing the degree did not necessarily lead to higher complexity, e.g. adding the degree $6$ term made the computation easier. A likely cause for this is that the dominating terms in the degree $6$ polynomial were simpler than those in the degree $5$ polynomial.