Real Root Isolation for Exp-Log Functions

Adam Strzebonski

Wolfram Research Inc.
100 Trade Centre Dr.
Champaign, IL 61820
U.S.A.

adams@wolfram.com
http://members.wolfram.com/adams

Abstract

We present a real root isolation procedure for univariate functions obtained by composition and rational operations from exp, log, and real constants. We discuss implementation of the procedure and give empirical results. The procedure requires the ability to determine signs of exp-log functions at simple roots of other exp-log functions. The currently known method to do this depends on Schanuel's conjecture.

Definition of exp-log functions

The set of exp-log functions is the smallest set of partial functions ℝ→ℝ containing exp, log, the identity function and the constant functions, closed under addition, multiplication and composition of functions.

Domains of exp-log functions

The domain D(f) of an exp-log function f is determined as follows:
• the domain of exp, the identity function and the constant functions is ℝ,
• the domain of log is (SubPlus[ℝ]),
• D(f+g)=D(f g)=D(f)⋂D(g),
• D(f(g))=(g^(-1))(D(f)).
In particular, D(f) is an open set and f is (C^∞) in D(f).

Examples of exp-log functions

The multiplicative inverse function:

inv:ℝ\{0}∋x→(1/x)=x exp (-log ((x^2)))∈ℝ

and the integer base radical functions:

(Subscript[rad, n]):(SubPlus[ℝ])∋x→(x^(1/n))=exp (((log(x))/n))∈ℝ, for n∈(SubPlus[ℤ])

are exp-log functions. Note, that the domain of (Subscript[rad, n]) does not include zero.

{x Exp[-Log[(x^2)]],Exp[((Log[x])/n)]}

Definition of roots

Let f:ℝ⊇D(f)→ℝ and let a∈D(f). a is a root of f iff 
• f(a)=0,
• (Subscript[∃, ϵ>0])(Subscript[∀, a-ϵ<x<a∨a<x<a+ϵ]) x∈D(f) ∧ x≠0 .

Finitness property of 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.

This property follows from the fact that the model-theoretic structure (Subscript[ℝ, exp])=(ℝ,<,0,1,+,·,exp) is o-minimal (Hovanskii 1980, Wilkie 1996).

The algorithm I will present computes domains and isolating intervals for real roots of exp-log functions. The proof of correctness of the algorithm does not use the o-minimality of (Subscript[ℝ, exp]) and hence it provides a new proof of the property.

Example

f:ℝ\{a}∋x→x-a-(Sqrt[(((x-a))^2)])∈ℝ

is an exp-log function which is zero at infinitely many points, but has no roots. Note, that the continuous extension of f to ℝ is not an exp-log function, since it is not differentiable at a.

Exp-log expressions

The set of exp-log expressions with coefficients in a computable field K⊆ℝ is defined recursively as follows:
• elements of K and the variable x are exp-log expressions,
• if f and g are exp-log expressions, so are f+g, f g and (f/g),
• if f is an exp-log expression, so are exp(f) and log(f).

Each exp-log expression represents an exp-log function, however the same function may be represented by many different expressions. From now on, when we refer to the domain, point values and roots of an exp-log expression, we will mean the domain, point values and roots of the corresponding exp-log function.

Two exp-log expressions are equivalent if they represent the same exp-log function. An exp-log function (expression) g extends an exp-log function (expression) f if D(f)⊆D(g) and (Subscript[∀, x∈D(f)]) f(x)=g(x).

Examples:
• exp(x)-exp(-x) is equivalent to ((exp(2x)-1)/(exp(x))).
• log((x^2)) extends 2 log(x).

For an exp-log expression f, f' denotes the expression obtained from f by applying the rules of differentiation with respect to x. The function represented by f' extends the derivative of the function represented by f.

Root constants

An exp-log root constant is an expression Root(f,a,b), where f is an exp-log expression, a<b are rational numbers, [a,b]⊆D(f), f (a) f (b)<0 and f' has a constant nonzero sign on (a,b). An exp-log root constant Root (f, a, b) represents the only root of f in (a,b). Form now on we will identify exp-log root expressions with the roots they represent.

Mathematica representation:

Reduce[E^x-3x-Log[x]==0,x,Reals]

r=%[[1,2]]

Precision[r]

N[r,100]

Root isolation algorithm

I will present a procedure which given an exp-log function finds its domain, roots and intervals on which the function is identically zero. The returned roots and interval endpoints are elements of K or exp-log root constants.

The assumption that K is a computable field implies that there is an algorithm for computation of signs of elements of K. However, the root isolation procedure requires the ability to determine signs of constants obtained by evaluation of exp-log functions at arbitrary exp-log root constants. At present, there is a known method to do this only when K is contained in the field of real elementary constants (Richardson, 1997). The proof that the method terminates relies on Schanuel' s conjecture. When K is not contained in the field of elementary constants, the root isolation procedure may fail if a call to the sign determination heuristic it uses fails.

Weak Fourier sequences

A sequence (Subscript[f, 1]),…,(Subscript[f, k]):I→ℝ of differentiable functions defined in an open interval I⊆ℝ is a weak Fourier sequence if 
• (Subscript[∀, 1≤i≤k-1])(Subscript[∀, x∈I]) sign((Subscript[f, i+1])(x))=sign((Subscript[f, i])'(x)),
• (Subscript[∀, x∈I]) (Subscript[f, k])(x)≠0.

Proposition: Let  (Subscript[f, 1]),…,(Subscript[f, k]):(a,b)→ℝ, where a,b∈ℝ⋃{-∞,∞}, be a weak Fourier sequence. Then for any 1≤i≤k
• (Subscript[f, i]) has at most k-i distinct roots in (a,b),
• (Subscript[∀, a≤x<b]) (Subscript[∃, x<c<b]) (Subscript[∀, x<y≤c]) (Subscript[f, i])(y)≠0,
• (Subscript[∀, a<x≤b]) (Subscript[∃, a<d<x]) (Subscript[∀, d≤y<x]) (Subscript[f, i])(y)≠0.

Notation: Let  (Subscript[f, 1]),…,(Subscript[f, k]):(a,b)→ℝ be a weak Fourier sequence. Then, for any 1≤i≤j≤k:
• sign((Subscript[f, i])((SuperPlus[x]))):=sign((Subscript[f, i])(c)), where a≤x<b and (Subscript[∀, x<y≤c]) (Subscript[f, i])(y)≠0,
• sign((Subscript[f, i])((SuperMinus[x]))):=sign((Subscript[f, i])(d)), where a<x≤b  and (Subscript[∀, d≤y<x]) (Subscript[f, i])(y)≠0,
• (Subscript[sgc, i,j])(x), for a<x<b, is number of sign changes in the sequence (Subscript[f, i])(x),…,(Subscript[f, j])(x) with terms equal to zero removed,
• (Subscript[sgc, i,j])((SuperPlus[x])), for a≤x<b, is number of sign changes in the sequence (Subscript[f, i])((SuperPlus[x])),…,(Subscript[f, j])((SuperPlus[x])),
• (Subscript[sgc, i,j])((SuperMinus[x])), for a<x≤b, is number of sign changes in the sequence (Subscript[f, i])((SuperMinus[x])),…,(Subscript[f, j])((SuperMinus[x])),
• sgc(x):=(Subscript[sgc, 1,k])(x), sgc((SuperPlus[x])):=(Subscript[sgc, 1,k])((SuperPlus[x])), sgc((SuperMinus[x])):=(Subscript[sgc, 1,k])((SuperMinus[x])).

Theorem: Let I⊆ℝ be an open interval and let (Subscript[f, 1]),…,(Subscript[f, k]):I→ℝ be a weak Fourier sequence. Then for any x∈I
• sgc((SuperPlus[x]))=sgc(x),
• sgc((SuperMinus[x]))=sgc(x)+r+2s, where r,s∈(SubPlus[ℤ])⋃{0},  (Subscript[f, 1])(x)=…=(Subscript[f, r])(x)=0 and (Subscript[f, r+1])(x)≠0.
Moreover, s=0 unless there is t>r+1 such that (Subscript[f, t])(x)=0.

Note that, unlike for polynomial derivative sequences, there is no relation between r and the multiplicity of the root of (Subscript[f, 1]) at x.

Semi Fourier sequences

A sequence ((Subscript[g, 1])/(Subscript[h, 1])),…,((Subscript[g, m])/(Subscript[h, m])):U→ℝ of functions defined in an open set U⊆ℝ is a semi Fourier sequence for f if U=D(f) and for all x∈U:
• (Subscript[h, i]) is continuous and (Subscript[h, i])(x)≠0, for 1≤i≤m,
• sign((((Subscript[g, 1])(x))/((Subscript[h, 1])(x))))=sign(f(x)),
• (Subscript[g, i]) is differentiable and sign((((Subscript[g, i+1])(x))/((Subscript[h, i+1])(x))))=sign((Subscript[g, i])'(x)), for 1≤i≤m-1,
• (Subscript[g, m]) is constant in each interval contained in U.

Corollary: Let  ((Subscript[g, 1])/(Subscript[h, 1])),…,((Subscript[g, m])/(Subscript[h, m])) be a semi Fourier sequence for f and let a<b be such that (a,b)⊆U. Functions (Subscript[h, 1]),…,(Subscript[h, m]) have constant nonzero signs (Subscript[s, 1]),…,(Subscript[s, m]) in (a,b). For 1≤i≤m put (Subscript[f, i]):=(Subscript[s, 1])·…·(Subscript[s, i])(Subscript[g, i]) so that sign((Subscript[f, i+1])(x))=sign((Subscript[f, i])'(x)) for all a<x<b. Let k be the largest index for which (Subscript[f, k]) is not identically zero on (a,b) or 0 if such index does not exist. If k>0, let 1≤p≤k. Then 
• (Subscript[f, k+1]), …,(Subscript[f, m]) are identically zero in (a,b),
• if k>0, (Subscript[f, k]) is nonzero and constant in (a,b),
• sign(f(x))=sign((Subscript[f, 1])(x)) for all x∈(a,b),
• (Subscript[sgc, p,k])((SuperPlus[a]))≥(Subscript[sgc, p,k])((SuperMinus[b])),
• if (Subscript[sgc, p,k])((SuperPlus[a]))=(Subscript[sgc, p,k])((SuperMinus[b])), then (Subscript[f, p]) has a constant nonzero sign in (a,b),
• if (Subscript[sgc, p,k])((SuperPlus[a]))=(Subscript[sgc, p,k])((SuperMinus[b]))+1, then (Subscript[f, p]) has exactly one simple root in (a,b).

Outline of the root isolation algorithm

Input: An exp-log expression f.
Output: {((Subscript[a, 1]),(Subscript[b, 1])),…,((Subscript[a, l]),(Subscript[b, l]))}, {((Subscript[c, 1]),(Subscript[d, 1])),…,((Subscript[c, m]),(Subscript[d, m]))}, {(Subscript[r, 1]),…,(Subscript[r, n])}, such that:
• (Subscript[a, i]), (Subscript[b, i]), (Subscript[c, j]), (Subscript[d, j]), (Subscript[r, k]) are elements of K, exp-log root constants, -∞, or ∞, 
• D(f)=((Subscript[a, 1]),(Subscript[b, 1]))⋃…⋃((Subscript[a, l]),(Subscript[b, l])),
• (f^(-1))(0)=((Subscript[c, 1]),(Subscript[d, 1]))⋃…⋃((Subscript[c, m]),(Subscript[d, m]))⋃{(Subscript[r, 1]),…,(Subscript[r, n])}.

(1) Find a semi Fourier sequence ((Subscript[g, 1])/(Subscript[h, 1])),…,((Subscript[g, t])/(Subscript[h, t])) for f, where (Subscript[g, i]) and (Subscript[h, i]) are exp-log expressions.
(2) Let H={(Subscript[h, 1]),…,(Subscript[h, t])}, let L be the set of all expressions that appear in f as arguments of log and let Q be the set of all expressions that appear in f as a denominator. Use the root isolation algorithm recursively to compute the domains and zero sets of H, L and Q.
(3) D(f)=((Subscript[a, 1]),(Subscript[b, 1]))⋃…⋃((Subscript[a, l]),(Subscript[b, l])) is the set where all elements of H and Q are defined and nonzero and all elements of L are defined and positive. 
(4) Use the corollary to isolate the roots of f in each ((Subscript[a, i]),(Subscript[b, i])).

Example

Let us follow the algorithm step-by-step.

Problem: Find the domain and the zero set of 
f(x)=((exp(x))/(1+x))-(x^2)log(x)

f=((Exp[x])/(1+x))-(x^2) Log[x]

Semi Fourier sequence

tmp=Together[f]

{(Subscript[g, 1]),(Subscript[h, 1])}={Numerator[tmp],Denominator[tmp]}

tmp=D[(Subscript[g, 1]),x]

{(Subscript[g, 2]),(Subscript[h, 2])}={Numerator[tmp],Denominator[tmp]}

tmp=D[(Subscript[g, 2]),x]

{(Subscript[g, 3]),(Subscript[h, 3])}={Numerator[tmp],Denominator[tmp]}

tmp=D[(Subscript[g, 3]),x]

{(Subscript[g, 4]),(Subscript[h, 4])}={Numerator[tmp],Denominator[tmp]}

tmp=D[(Subscript[g, 4]),x]

The log term is gone. Put the terms over a common denominator.

tmp=Together[tmp]

{(Subscript[g, 5]),(Subscript[h, 5])}={Numerator[tmp],Denominator[tmp]}

tmp=D[(Subscript[g, 5]),x]

{(Subscript[g, 6]),(Subscript[h, 6])}={Numerator[tmp],Denominator[tmp]}

tmp=D[(Subscript[g, 5]),x]

{(Subscript[g, 6]),(Subscript[h, 6])}={Numerator[tmp],Denominator[tmp]}

tmp=D[(Subscript[g, 6]),x]

Divide out by a term that is always positive.

tmp=Cancel[tmp/E^x]

{(Subscript[g, 7]),(Subscript[h, 7])}={Numerator[tmp],Denominator[tmp]}

tmp=D[(Subscript[g, 7]),x]

{(Subscript[g, 8]),(Subscript[h, 8])}={Numerator[tmp],Denominator[tmp]}

tmp=D[(Subscript[g, 8]),x]

{(Subscript[g, 9]),(Subscript[h, 9])}={Numerator[tmp],Denominator[tmp]}

The domain

H=((Subscript[h, #])&)/@Range[9]

L=First/@Cases[f,_Log,{0,∞}]

Q=First/@Cases[f,Power[_,_?Negative],{0,∞}]

The domain (in general a recursive call to the exp-log root isolation algorithm).

Reduce[And@@Thread[Join[H,Q]≠0]&&And@@Thread[L>0],x,Reals]

Root isolation in D(f)=(0,∞)

Signs of denominators.

Sign[H/.x->1]

All denominators are positive, hence (Subscript[f, i])=(Subscript[g, i]).

F=((Subscript[g, #])&)/@Range[9]

sign((Subscript[f, i])((SuperPlus[0]))):

(Subscript[s, 1])=Sign[Limit[F,x->0]]

sign((Subscript[f, i])((SuperMinus[∞]))):

(Subscript[s, 2])=Sign[Limit[F,x->∞]]

The smallest p for which (Subscript[sgc, p,9])((SuperPlus[0]))=(Subscript[sgc, p,9])((SuperMinus[∞]))+1 is 6. Find the root of (Subscript[f, 6]).

Reduce[F[[6]]==0 &&x>0,x,Reals]

(Subscript[r, 1])=%[[2]]

(Subscript[f, 6]) is nonzero in (0,(Subscript[r, 1])) and in ((Subscript[r, 1]),∞).

Interval (0, (Subscript[r, 1]))

sign((Subscript[f, i])((SuperMinus[(Subscript[r, 1])]))) for 1≤i≤6:

(Subscript[s, 3])=Sign[Append[Take[F,5]/.x->(Subscript[r, 1]),F[[6]]/.x->(Subscript[r, 1])/2]]

Take[(Subscript[s, 1]),6]

The smallest p for which (Subscript[sgc, p,6])((SuperPlus[0]))=(Subscript[sgc, p,6])((SuperMinus[(Subscript[r, 1])]))+1 is 5. Find the root of (Subscript[f, 5]).

(Subscript[r, 2])=Reduce[F[[5]]==0 &&0<x<(Subscript[r, 1]),x,Reals][[2]]

(Subscript[f, 5]) is nonzero in (0,(Subscript[r, 2])) and in ((Subscript[r, 2]),(Subscript[r, 1])).

Interval (0, (Subscript[r, 2]))

sign((Subscript[f, i])((SuperMinus[(Subscript[r, 2])]))) for 1≤i≤5:

Sign[Append[Take[F,4]/.x->(Subscript[r, 2]),F[[5]]/.x->(Subscript[r, 2])/2]]

Take[(Subscript[s, 1]),5]

(Subscript[sgc, 1,5])((SuperPlus[0]))=(Subscript[sgc, 1,5])((SuperMinus[(Subscript[r, 2])]))=2, hence there are no roots of f in (0,(Subscript[r, 2])).

Interval ((Subscript[r, 2]),(Subscript[r, 1]))

sign((Subscript[f, i])((SuperPlus[(Subscript[r, 2])]))) for 1≤i≤5:

Sign[Append[Take[F,4]/.x->(Subscript[r, 2]),F[[5]]/.x->((Subscript[r, 2])+(Subscript[r, 1]))/2]]

Take[(Subscript[s, 3]),5]

(Subscript[sgc, 1,5])((SuperPlus[(Subscript[r, 2])]))=(Subscript[sgc, 1,5])((SuperMinus[(Subscript[r, 1])]))=1, hence there are no roots of f in ((Subscript[r, 2]),(Subscript[r, 1])).

Interval ((Subscript[r, 1]),∞)

sign((Subscript[f, i])((SuperPlus[(Subscript[r, 1])]))) for 1≤i≤6:

(Subscript[s, 4])=Sign[Append[Take[F,5]/.x->(Subscript[r, 1]),F[[6]]/.x->2(Subscript[r, 1])]]

Take[(Subscript[s, 2]),6]

The smallest p for which (Subscript[sgc, p,6])((SuperPlus[(Subscript[r, 1])]))=(Subscript[sgc, p,6])((SuperMinus[∞]))+1 is 3. Find the root of (Subscript[f, 3]).

(Subscript[r, 3])=Reduce[F[[3]]==0 &&x>(Subscript[r, 1]),x,Reals][[2]]

(Subscript[f, 3]) is nonzero in ((Subscript[r, 1]),(Subscript[r, 3])) and in ((Subscript[r, 3]),∞).

Interval ((Subscript[r, 1]),(Subscript[r, 3]))

sign((Subscript[f, i])((SuperMinus[(Subscript[r, 3])]))) for 1≤i≤3:

Sign[Append[Take[F,2]/.x->(Subscript[r, 3]),F[[3]]/.x->((Subscript[r, 1])+(Subscript[r, 3]))/2]]

Take[(Subscript[s, 4]),3]

(Subscript[sgc, 1,3])((SuperPlus[(Subscript[r, 1])]))=(Subscript[sgc, 1,3])((SuperMinus[(Subscript[r, 3])]))+1, hence there is exactly one root of f in ((Subscript[r, 1]),(Subscript[r, 3])).

(Subscript[ans, 1])=Reduce[f==0 &&(Subscript[r, 1])<x<(Subscript[r, 3]),x,Reals][[2]]

Interval ((Subscript[r, 3]),∞)

sign((Subscript[f, i])((SuperPlus[(Subscript[r, 3])]))) for 1≤i≤3:

Sign[Append[Take[F,2]/.x->(Subscript[r, 3]),F[[3]]/.x->2(Subscript[r, 3])]]

Take[(Subscript[s, 2]),3]

(Subscript[sgc, 1,3])((SuperPlus[(Subscript[r, 3])]))=(Subscript[sgc, 1,3])((SuperMinus[∞]))+1, hence there is exactly one root of f in ((Subscript[r, 3]),∞).

(Subscript[ans, 2])=Reduce[f==0 &&x>(Subscript[r, 3]),x,Reals][[2]]

Answer

{(Subscript[ans, 1]),(Subscript[ans, 2])}

Examples (Gruntz)

Reduce[(E^x) ((E^(-(E^(-x))+(1/x)))-(E^(1/x)))+5==0,x,Reals]

Reduce[(E^(-(E^(E^x))+(E^(E^((E^(-x))+x)))))-(10^67)==0,x,Reals]

Reduce[((((3^x)+(5^x)))^(1/x))-10==0,x,Reals]

Reduce[(E^(-x+(x/((E^x) ((E^(-x))+(E^(-((2 (x^2))/(1+x)))))))))-(7/10)==0,x,Reals]

Examples (MathWorld)

Smarandache constant.

Reduce[(127^x)-(113^x)-1==0,x,Reals]

Laplace limit constant.

Reduce[((x Exp[(Sqrt[1+(x^2)])])/(1+(Sqrt[1+(x^2)])))-1==0,x,Reals]

Fixed point of ArcCoth

fps=Reduce[ArcCoth[x]==x,x,Reals]

Show[{Plot[{ArcCoth[x],x},{x,-3,3}],Graphics[{Red,PointSize[Medium],Point[{x,x}]/.{ToRules[fps]}}]}]

Inflection points of Einstein's (Subscript[E, 1]) function

Einstein's (Subscript[E, 1]) function.

(Subscript[E, 1])=(((x^2) Exp[x])/(((Exp[x]-1))^2));

infl=Reduce[(∂_({x,2}))(Subscript[E, 1])==0,x,Reals]

Show[{Plot[(Subscript[E, 1]),{x,-5,5}],Graphics[{Red,PointSize[Medium],Point[{x,(Subscript[E, 1])}]/.{ToRules[infl]}}]}]

Intersection of Einstein's (Subscript[E, 4]) function with hyperbola

Einstein's (Subscript[E, 4]) function.

(Subscript[E, 4])=(x/(Exp[x]-1))-Log[1-Exp[-x]];

intp=Reduce[(Subscript[E, 4])==(1/x),x,Reals]

Show[{Plot[{(Subscript[E, 4]),1/x},{x,0,3}],Graphics[{Red,PointSize[Medium],Point[{x,(Subscript[E, 4])}]/.{ToRules[intp]}}]}]

Planck's radiation function

P=(15/((π^4) ((x^5) (Exp[(1/x)]-1))));

Maximum of P.

max=Reduce[(∂_x)P==0,x,Reals]

Inflection points of P.

infl=Reduce[(∂_({x,2}))P==0,x,Reals]

Show[{Plot[P,{x,0,1}],Graphics[{PointSize[Medium],{Red,Point[{x,P}]/.{ToRules[max]}},{Green,Point[{x,P}]/.{ToRules[infl]}}}]}]

High degree sparse polynomial

Reduce[(x^1000000)-2 (x^777777)+3 (x^123456)-4 (x^54321)+5 (x^6789)-6 (x^555)+7 (x^100)-8 (x^21)+9 x-10==0,x,Reals]

Algebraic function involving high degree radicals

Reduce[2 (x^(123451/67890))-(x^2)+4 (Sqrt[x])-4 x-(9/8)==0,x,Reals]

Function involving non-rational real powers

Reduce[(x^π)-(x^(x^(Sqrt[2])))-(Sqrt[3]) x+(2^(1/3))==0,x,Reals]

Function with a quadruple root at 1 and two simple roots

Reduce[-1+(E^(-1+x))-(4/5) (-1+x)-(8/25) (((-1+x))^2)-(51/250) (((-1+x))^3)+Log[5]-Log[3+(E^(((-1+x))^2))+x]==0,x,Reals]

Examples with many roots

{(Subscript[F, 1]),(Subscript[F, 2]),(Subscript[F, 3])}={Exp[ChebyshevT[7,x]]-Log[ChebyshevU[7,x]+10],Exp[ChebyshevT[10,x]]-Log[ChebyshevU[10,x]+10],Exp[ChebyshevT[10,x]]-9 Log[ChebyshevU[10,x]+10]-(x^2)+(389/20)};

Plot[Evaluate[{(Subscript[F, 1]),(Subscript[F, 2]),(Subscript[F, 3])}],{x,-1.1,1.1}]

Length[Reduce[#==0,x,Reals]]&/@{(Subscript[F, 1]),(Subscript[F, 2]),(Subscript[F, 3])}