Approximating the pth Root by Composite Rational Functions
Evan S. Gawlik Yuji Nakatsukasa

Concentration tuned tetragonal strain in alloys: application to magnetic anisotropy of FeNi_1-xCo_x
    Vladimir P. Antropov
    December 30, 2023
===================================================================================================




A landmark result from rational approximation theory states that x^1/p on [0,1] can be approximated by a type-(n,n) rational function with root-exponential accuracy. 
Motivated by the recursive optimality property of Zolotarev functions (for the square root and sign functions), we investigate approximating x^1/p by composite rational functions of the form 
r_k(x, r_k-1(x, r_k-2( ⋯ (x,r_1(x,1)) ))). 
While this class of rational functions ceases to contain the minimax (best) approximant for p≥ 3, we show that it achieves approximately pth-root exponential convergence with respect to the degree. 
Moreover, crucially, the convergence is doubly exponential with respect to the number of degrees of freedom, suggesting that composite rational functions are able to approximate x^1/p and related functions (such as |x| and the sector function) with exceptional efficiency. 












§ INTRODUCTION


Composing rational functions is an efficient way of generating a rational function r(x)=r_k(⋯ r_2(r_1(x))) of high degree: if each r_i is of type (m,m), then r is of type (m^k,m^k). By choosing each r_i appropriately, one can often obtain a function r that approximates a desired function in a wide domain of interest. 

There is no reason to expect—and it is generally not true—that 
r_k(⋯ r_2(r_1(x))) can express the minimax rational approximant of a given type, say (m^k,m^k), to a given function. 
However, building upon Rutishauser <cit.> and Ninomiya <cit.>, Nakatsukasa and Freund <cit.> show a remarkable property of the best rational approximants to the function sign(x)=x/|x| on [-1,-δ]∪ [δ,1] for 0<δ<1 (called Zolotarev functions): appropriately composing Zolotarev functions gives another Zolotarev function of higher degree. In other words, the class of composite rational functions r(x)=r_k(⋯ r_2(r_1(x))), with each r_i of type (m,m), contains the type-(m^k,m^k) minimax approximant to the sign function.
Moreover, for a fixed δ, the convergence of Zolotarev functions is exponential in the degree. Since the degree is m^k, and the number of parameters necessary to express r is d≈ 2km, 
it follows that the convergence is exp(-m^k) = exp(-exp(C d)), a double-exponential convergence rate. This is so powerful that choosing m=17 and k=2 (one composition, i.e., two iterations) is enough to obtain convergence to machine precision in double precision arithmetic, with error below 10^-15. 

Functions related to the sign function, such as 
|x| (via |x|=x/(x)) and √(x) 
(via |x|≈ p(x^2)/q(x^2) then √(x)≈ p(x)/q(x)) can similarly be approximated by composite rational functions. Gawlik <cit.> does this for the square root and shows that 
a composite rational function yields the minimax rational approximant (in the relative sense) on intervals [δ,1] ⊂ (0,1], and that the approximation extends far into the complex plane.  This observation generalizes earlier work on rational approximation of the square root with optimally scaled Newton iterations <cit.>.
Moreover, an extension was derived in <cit.>, which shows that the pth root can be approximated efficiently on intervals [δ,1] ⊂ (0,1], although not with minimax quality. 

Clearly, in the above papers the origin is excluded from the domain, as the functions have a singularity at x=0. However, a landmark result from rational approximation theory <cit.> states that the best rational approximant (in the absolute sense) of x^β (for any real β >0) on [0,1] can be approximated by a type-(n,n) rational function with root-exponential accuracy. One might wonder, can this be done with a composite rational function? This is the question we address in this paper.  We focus on the case in which β = 1/p with p ≥ 2 an integer.

We show that a rational function of the form r(x) = r_k(x, r_k-1(x, r_k-2( ⋯ (x,r_1(x,1)) ))) can approximate x^1/p on [0,1] with superalgebraic accuracy, with close to pth root-exponential convergence. Moreover—and crucially—the convergence is doubly exponential with respect to the number of degrees of freedom.  That is, the error is O(exp(-c_1exp(c_2d))) for some constants c_1,c_2>0, where d is the number of parameters needed to express the rational function.  By “number of parameters” we mean d=∑_i=1^k m_i+ℓ_i+1 if r_i has type (m_i,ℓ_i) for i=1,2,…,k, so that d reflects the cost of evaluating r at a matrix argument. 

Clearly, our result implies that any rational power of x can be 
approximated by a composite rational function. Moreover, 
since |r(x)-x^1/p|≤ϵ on [0,1] implies 
|r(x/s)- (x/s)^1/p|≤ϵ on [0,s] for any s>0, hence 
|s^1/pr(x/s)- x^1/p|≤ s^1/pϵ, 
our results also show that any rational power can be approximated efficiently on [0,s] by a composite rational function. In addition, our approximants to x^1/p immediately lead to approximants to the p-sector function sect_p(z)=z/(z^p)^1/p.


More generally, we think composite (rational) functions are a powerful tool in approximation theory, and we regard this as a contribution towards demonstrating their effectiveness and practicality. 
Indeed, one might say they are already used extensively in scientific computing: 


  * Composite rational functions are implicitly employed in most algorithms for computing matrix functions <cit.>, in which approximating a function on the spectrum of the matrix is required. For the pth root, a standard algorithm <cit.> employs Newton's method, which ultimately approximates A^1/p with a sequence of rational functions f_k of A given recursively by f_k+1(x) = 1/p((p-1)f_k(x)+x/f_k(x)^p-1), f_0(x)=1.   The function f_k is composite rational and similar to the approximants we use, but not the same (it is unscaled), and it exhibits exponential rather than double-exponential convergence on [0,1].  
Generally speaking, Newton's method for computing a matrix function f(A) (or more generally for various nonlinear problems, e.g. rootfinding) can often be interpreted as approximating f(A) (or the solution) by a composite rational function of A.

  * The rapidly growing subject of deep learning is based on composing a large number of nonlinear activation functions <cit.>. 




  
Summary of Results.
To summarize our results, let us introduce some terminology.
We say that a univariate rational function r(x)=p(x)/q(x) is of type (m,ℓ) if p and q are polynomials of degrees at most m and ℓ, respectively.
We denote the set of all such rational functions by ℛ_m,ℓ.
We say that a bivariate rational function r(x,y) is of type (m,ℓ) if r(x,x) is of type (m,ℓ).  
We say that a univariate rational function r is (k,m,ℓ)-composite if r is a composition of k rational functions r_i(x,y), i=1,2,…,k, each of type (m,ℓ):

    r(x) = r_k(x, r_k-1(x, r_k-2( ⋯ (x,r_1(x,1))))).


Here is the main result of this paper. 

Let p ≥ 2 be an integer.  There exists a positive constant N depending on p such that for every integer n ≥ N, there exists a (⌊log_p n ⌋ +1,p,p-1)-composite rational function r of type (n,n-1) such that

    max_x ∈ [0,1] |r(x)-x^1/p| ≤exp( -b n^c ),

where b>0 is a constant depending on p and

    c = log(p/p-1) log 2 /log(2p/p-1) log p .


Note that when p=2, c=1/2, and as p →∞, c ∼1/plog p.

Let us comment on the theorem. The bound (<ref>) shows that by using a (⌊log_p n ⌋ +1,p,p-1)-composite rational function we can approximate the pth root with “1/cth root”-(nearly pth root) exponential accuracy with respect to the degree, which is suboptimal unless p=2 (in which case a composite rational function on [δ,1] is optimal in the relative sense). 

However, the result is still striking in the following sense: the number of degrees of freedom used to express r is just O(pk) for n≈ p^k (see below (<ref>)), and therefore with respect to the degrees of freedom d, 
the convergence is 

    max_x ∈ [0,1] |r(x)-x^1/p| ≤exp( -b p^c̃d )  ,
 indicating a double-exponential convergence with respect to d. 

As a byproduct of our analysis, we will obtain analogous results for composite rational approximation of the p-sector function sect_p(z)=z/(z^p)^1/p on the set S_p ⊂ℂ given by

    S_p = { x e^2π i j/p| x ∈ [0,1],   j ∈{1,2,…,p}}.

We will also consider the subset S_p,α of S_p excluding the origin 

    S_p,α = { x e^2π i j/p| x ∈ [α,1],   j ∈{1,2,…,p}}.


We say that a (k,m,ℓ)-composite rational function (<ref>) is pure if the functions r_j(x,y) appearing in (<ref>) are univariate:

    r(x) = r_k(r_k-1( ⋯ (r_1(x)))).

 
Let p ≥ 2 be an integer, and α∈ (0,1).  
There exists a positive constant N depending on p such that for every integer n ≥ N, there exist pure (⌊log_p n ⌋,1,p)-composite rational functions r and q of type (n-p+1,n) such that

    max_z ∈ S_p|z ( r(z) - sect_p(z) )| ≤exp( -b n^c ),

where b and c are as in Theorem <ref>, and 

    max_z ∈ S_p,α| q(z) - sect_p(z) | ≤exp( -b n^c ),

where b>0 depends on α and p, and c = log 2/log p.  

It is worth noting that the two rational functions r, q are 
generally different—they coincide for a particular value of α.
The error in (<ref>) is measured in a weighted norm, which is natural in view of the fact that sect_p(z) is discontinuous at z=0. When p=2 and z ∈ S_2, zsect_p(z)=|z| and 
c=1/2, so (<ref>) recovers the root-exponential convergence of rational approximants to |x| on [-1,1] <cit.>. 
By contrast, (<ref>) shows that a better bound holds for the absolute error if one excludes the neighborhood of the origin. When p=2, c=1 and (<ref>) recovers the exponential convergence of Zolotarev functions to the sign function on [-1,-α]∪ [α,1] <cit.>. Our analysis will show that b decays like a negative power of log1/α as α→ 0.



  
Organization.  
This paper is organized as follows.  In Section <ref>, we review some theory from <cit.> concerning composite rational approximants of the pth root on positive real intervals.  In Section <ref>, we study the behavior of these approximants near the origin.  We then prove Theorems <ref> and <ref> in Section <ref>, and we illustrate our results numerically in Section <ref>.



§ COMPOSITE RATIONAL APPROXIMATION OF THE PTH ROOT
 
To approximate x^1/p on an interval [α^p,1] ⊂ (0,1], Gawlik <cit.> considers the recursively defined rational function

    f_k+1(x)    = f_k(x) _m,ℓ(x/f_k(x)^p,α_k, · ),       f_0(x)   =1, 
    α_k+1   = α_k /_m,ℓ(α_k^p,α_k, · ) ,      α_0    = α,

where _m,ℓ(x,α, · ) is (a rescaling of) the relative minimax rational approximant of type (m,ℓ) ∈ℕ_0 ×ℕ_0 ∖{(0,0)} on the interval [α^p,1]:

    r̂_m,ℓ(x,α, · ) = ( 1+α/2α) r_m,ℓ(x,α, · ),

where

    r_m,ℓ( · , α,  · ) = _r ∈ℛ_m,ℓmax_x ∈ [α^p,1]| r(x)-x^1/p/x^1/p|.

Gawlik shows that f_k(x) is a rapidly convergent approximant to the pth root on [α^p,1]. With k recursions, the maximum relative error |f_k(x)-x^1/p|/|x^1/p| on [α^p,1] decays double exponentially in k: it is bounded above by c_1 exp(-c_2(m+ℓ+1)^k) for some c_1,c_2>0 depending on m, ℓ, p, and α.  Importantly, these constants depend very weakly on α; the analysis below will implicitly show that when (m,ℓ)=(1,0), c_1 is independent of α and c_2 decays like a negative power of log1/α as α→ 0, just like b in (<ref>).

Given that (<ref>) is an approximant on [α^p,1], which is an interval that excludes the singularity at x=0, 
a natural question arises: can we approximate on [0,1]? 
Intuitively, 
the function is still continuous at x=0 (unlike e.g. the sign or sector function) with 0^1/p=0, and hence it is possible to approximate x^1/p on the whole interval [0,1]. Indeed Stahl <cit.> shows that x^1/p on [0,1] can be approximated by a type-(n,n) rational function with root-exponential accuracy (we refer to <cit.> for general results on classical rational approximation theory). 
Can a highly efficient rational approximant be constructed based on recursion as in (<ref>)? 
It is important to note that we will necessarily switch to the (more natural) metric of absolute error |r(x)-x^1/p| rather than the relative error |r(x)-x^1/p|/|x^1/p| for this purpose. 

It turns out that the rational function (<ref>) does a good job approximating on [0,1], when α is chosen carefully: when it is too small, the error is large on [α^p,1] (in fact it is maximal at x=1 <cit.>).
Conversely if α is too large, the error is large on [0,α^p] (in fact it is O(α) at x=0, as we show below). 
A major task undertaken in what follows is to choose  α so that the convergence is optimized, in that the error on [0,α^p] and [α^p,1] are balanced to be approximately the same. 

Our analysis will focus on the lowest-order version of the iteration (<ref>-<ref>),  obtained by choosing (m,ℓ)=(1,0).  It is shown in <cit.> (and elsewhere <cit.>) that for this choice of m and ℓ, 

    _1,0(x,α, · )=1/p((p-1)μ(α)+x/μ(α)^p-1),    μ(α) =(α-α^p/(p-1)(1-α))^1/p.

Thus, when (m,ℓ)=(1,0), the iteration (<ref>-<ref>) reads

    f_k+1(x)    = 1/p( (p-1) μ(α_k) f_k(x) + x/μ(α_k)^p-1 f_k(x)^p-1),       f_0(x)    = 1, 
    α_k+1   = p α_k/(p-1)μ(α_k) + μ(α_k)^1-pα_k^p ,      α_0    = α.

Note that f_k is (k,p,p-1)-composite since it is of the form (<ref>) with

    r_j(x,y) = 1/p(  (p-1) μ(α_j-1)^p y^p + x /μ(α_j-1)^p-1 y^p-1)

for each j.  It follows from this observation and an inductive argument that f_k has type (p^k-1,p^k-1-1) for each k ≥ 1.


We rely heavily on this explicit expression for the particular case (m,ℓ)=(1,0), as it lets us analyze the functions in detail, which leads to a constructive proof for Theorem <ref>. 
We note that using larger values of (m,ℓ) may result in faster convergence, in particular a larger exponent c than (<ref>). In view of (<ref>), the convergence is still doubly exponential, with an improved constant c̃. However, we do not expect the improvement would be significant.

Moreover, composing low-degree rational functions is an extremely efficient way to construct high-degree rational functions of matrices, and we suspect that our choice (m,ℓ)=(1,0) would give the fastest convergence in terms of the number of matrix operations needed to evaluate r at a matrix argument.



§ BOUNDING THE ERROR ON [0,Α^P]
 

In this section, we analyze the absolute error committed by the function f_k defined by (<ref>)–(<ref>) on the interval [0,α^p].  It will be convenient to consider not f_k but the scaled function 

    f_k(x) = 2α_k/1+α_k f_k(x),

which has the property that <cit.>

    max_x ∈ [α^p,1]f_k(x) - x^1/p/x^1/p = -min_x ∈ [α^p,1]f_k(x) - x^1/p/x^1/p  = 1-α_k/1+α_k∈ (0,1).



We will prove the following estimate.

 
Let α∈ (0,1).  The function f_k defined by (<ref>)–(<ref>) and (<ref>) satisfies

    max_x ∈ [0,α^p] | f_k(x) - x^1/p  | ≤ 2α

for every k ≥ 0.


Experiments suggest that the bound (<ref>) could be improved to <α for k large enough, but this does not affect what follows in any significant way.

We will prove Theorem <ref> by a series of lemmas.  Let

    g_k(x) = x/f_k(x^p).

Note that g_0(x) = x and

    g_k+1(x) = x/ f_k(x^p) r̂_1,0( x^p/f_k(x^p)^p, α_k,  · )  = g_k(x)/r̂_1,0(g_k(x)^p, α_k, · )  = ŝ(g_k(x), α_k),

where

    ŝ(x,α) = x/r̂_1,0(x^p,α, · ) = px/(p-1)μ(α) + μ(α)^1-p x^p.

Also let

    H(α) = ŝ(α,α) = p α/(p-1)μ(α) + μ(α)^1-pα^p,

so that α_k+1 = H(α_k).

 
For every α∈ (0,1) and every x ∈ [0,α],

    0 ≤ x ŝ'(x,α) ≤ŝ(x,α) ≤ H(α),

where the prime denotes differentiation with respect to x.


A short calculation shows that

    x ŝ'(x,α) = w(x) ŝ(x,α),

where 

    w(x) =  (p-1)(1 - (x/μ(α))^p ) / (p-1) + (x/μ(α))^p  .

Since 0 ≤ w(x) ≤ 1 for every x ∈ [0,μ(α)], it follows that 

    0 ≤ x ŝ'(x,α) ≤ŝ(x,α),    x ∈ [0,μ(α)].

In particular, the above inequalities hold on [0,α] ⊂ [0,μ(α)], and ŝ(x,α) is nondecreasing on [0,α].  Thus,

    ŝ(x,α) ≤ŝ(α,α) = H(α),    x ∈ [0,α].



Now let α∈ (0,1) be fixed.
 
For every x ∈ [0,α] and every k ≥ 0,

    0 ≤ x g_k'(x) ≤ g_k(x) ≤α_k.



Since g_0(x)=x and α_0=α, the above inequalities hold when k=0.  Assume that they hold for some k ≥ 0.
Observe that

    x g_k+1'(x) = x g_k'(x) ŝ'(g_k(x), α_k) .

Since g_k(x) ∈ [0,α_k] for x ∈ [0,α], Lemma <ref> implies that ŝ'(g_k(x), α_k) ≥ 0.  It follows from this and our inductive hypothesis that xg_k+1'(x) ≥ 0 for x ∈ [0,α].  In addition, since x g_k'(x) ≤ g_k(x) and g_k(x) ŝ'(g_k(x),α_k) ≤ŝ(g_k(x),α_k),

    x g_k+1'(x) ≤ŝ(g_k(x),α_k) = g_k+1(x).

Finally, since ŝ(g_k(x),α_k) ≤ H(α_k) = α_k+1, it follows that g_k+1(x) ≤α_k+1.



For every x ∈ [0,α^p] and every k ≥ 0,

    0 < f_k(x) ≤α(1+ε_k),   ε_k = 1-α_k/1+α_k.



We first note that f_k is positive and nondecreasing on [0,α^p].  Indeed, differentiating the relation

    f_k(x^p) = x/g_k(x)

gives

    p x^p-1 f_k'(x^p) = g_k(x) - x g_k'(x)/g_k(x)^2,

so Lemma <ref> implies that f_k'(x^p) ≥ 0 for every x ∈ [0,α].  
Evaluating the recursion (<ref>) at x=0 gives

    f_k+1(0) = f_k(0) ( p-1/p) μ(α_k),    f_0(0) = 1,

so f_k(0)>0 for every k.
Since f_k(x) is a positive multiple of f_k(x), it follows that 0 < f_k(x) ≤f_k(α^p) for every x ∈ [0,α^p].  Finally, taking x=α^p in (<ref>) gives f_k(α^p) ≤α(1+ε_k).


By the lemma above,

    |f_k(x) - x^1/p| ≤max{ |f_k(x)|, |x^1/p| }≤max{α(1+ε_k),  α} = α(1+ε_k) ≤ 2 α,    x ∈ [0,α^p],

so

    max_x ∈ [0,α^p] |f_k(x) - x^1/p| ≤ 2α.

This completes the proof of Theorem <ref>.

An estimate for the absolute error on [0,1] is now immediate: Combining the above theorem, (<ref>), and the fact that x^1/p≤ 1 for x ∈ [0,1], we see that

    max_x ∈ [0,1] | f_k(x) - x^1/p  | ≤max{ 2α, 1-α_k/1+α_k}.




 §.§ Sector function approximation

We note that the function g_k in (<ref>)
approximates the p-sector function _p(z)=z/(z^p)^1/p (this observation appeared in <cit.>), and 
g_k is a pure composite rational function of the form 
g_k(z) = r_k(r_k-1(⋯ r_2(r_1(z)))). 
In fact it is (k,1,p)-composite, and an inductive argument shows that it has type (p^k-p+1,p^k).
In the p=2 case, this reduces to Zolotarev's best rational 
approximant to the sign function of type (2^k-1,2^k).  That is, as in the square root approximation, the minimax rational approximant is contained in the class of (here purely) composite rational functions. 

Below we derive estimates for the maximum weighted error |z(g_k(z) - sect_p(z))| on the sets S_p,S_p,α⊂ℂ defined in (<ref>) and (<ref>).  As before, it will be convenient to work not with g_k(z) but with the rescaled function

    g_k(z) = 2/1+α_k g_k(z) =  4α_k/(1+α_k)^2z/f_k(z^p).

As shown in <cit.>,  the relative error g_k(z)-sect_p(z)/sect_p(z) is real-valued and equioscillates on each line segment {z∈ℂ| e^-2π ij/pz∈[α,1] }, j=0,1,…,p-1. 
Note that here the relative and absolute errors are the same in modulus. The asymptotic convergence rate on S_p,α was analyzed in <cit.>. Here we quantify the non-asymptotic convergence on S_p.

For every k ≥ 0,

    max_z ∈ S_p |z (g_k(z) - sect_p(z))| ≤max{α, 1-α_k/1+α_k},

and 

    max_z ∈ S_p,α |g_k(z) - sect_p(z)| ≤1-α_k/1+α_k .



Let z = x^1/p e^2π i j/p with x ∈ [0,1] and j ∈{1,…,p}.  Since g_k(z) = e^2π i j/pg_k(x^1/p) and sect_p(z) = e^2π i j/p, we have

    |z (g_k(z) - sect_p(z))| = |x^1/p (g_k(x^1/p) - 1)|.

If x ∈ [0,α^p], then Lemma <ref> implies that 0 ≤g_k(x^1/p) ≤2α_k/1+α_k < 1, so

    |x^1/p (g_k(x^1/p) - 1)| ≤ x^1/p≤α,    x ∈ [0,α^p].

On the other hand, if x ∈ [α^p,1], then

    |x^1/p (g_k(x^1/p) - 1)| ≤ |g_k(x^1/p) - 1| = | 4α_k/(1+α_k)^2x^1/p/f_k(x) - 1|.

By (<ref>),

    f_k(x)/x^1/p∈[ 1 - ( 1-α_k/1+α_k), 1 + ( 1-α_k/1+α_k) ] = [ 2α_k/1+α_k, 2/1+α_k] ,    x ∈ [α^p,1],

so

    x^1/p/f_k(x)∈[ 1+α_k/2, 1+α_k/2α_k] ,    x ∈ [α^p,1],

and hence

    4α_k/(1+α_k)^2x^1/p/f_k(x) - 1 ∈[ -1-α_k/1+α_k, 1-α_k/1+α_k] ,    x ∈ [α^p,1].

It follows that 

    |x^1/p (g_k(x^1/p) - 1)| ≤1-α_k/1+α_k,    x ∈ [α^p,1].

For (<ref>), we simply start from the second expression in (<ref>) and use (<ref>).




§ PROOF OF THEOREMS <REF> AND <REF>
 
To examine the convergence of the recursion (<ref>)-(<ref>) on [0,1], we first ask the question: given ϵ>0, what values of k and α are needed to get an error ϵ?  In view of (<ref>), we must choose α≤ϵ/2 and k large enough so that 1-α_k/1+α_k≤ϵ.

To determine k, we select a constant α^* ∈ (1/e,1) (depending on p) and split the convergence of α_k→ 1 into three stages: 


  * Find k_1 such that α_k_1≥1/e. 

  * Find k_2 such that α_k_1+k_2≥α^*. 

  * Find k_3 such that 1-α_k_1+k_2+k_3/1+α_k_1+k_2+k_3≤ϵ. 

Clearly, the second stage is independent of ϵ and α_0, so k_2 is a constant (depending on p). 

Our choice of α^* is described in the following lemma.

 
There exists a constant α^* ∈ (0,1) depending on p such that  

    1-H(α)/1+H(α)≤p/2(1-α/1+α)^2

for every α∈ [α^*,1].


It is proven in <cit.> that the iteration (<ref>) generates an increasing sequence {α_k}_k=0^∞ satisfying lim_k →∞α_k = 1 and

    1-α_k+1/1+α_k+1 = C(m,ℓ,p) ( 1-α_k/1+α_k)^m+ℓ+1 + o ( ( 1-α_k/1+α_k)^m+ℓ+1),

where

    C(m,ℓ,p) = p^m+ℓ+1 m! ℓ! (1/p)_ℓ+1 (1-1/p)_m / 2^m+ℓ (m+ℓ+1)!(m+ℓ)!

and (β)_m denotes the rising factorial (the Pochhammer symbol): (β)_m = β(β+1)(β+2)⋯(β+m-1).  Since C(1,0,p) = p-1/4, this implies that the iteration (<ref>) with (m,ℓ)=(1,0) (i.e., the iteration (<ref>)) generates {α_k}_k=0^∞ satisfying

    1-α_k+1/1+α_k+1 = ( p-1/4) ( 1-α_k/1+α_k)^2 + o ( ( 1-α_k/1+α_k)^2 ).

In other words,

    . 1-H(α)/1+H(α)/ (1-α/1+α)^2 .  →p-1/4,  as α↑ 1.

It follows that the above ratio is bounded by p/2 for α close enough to 1.

Without loss of generality, we assume 

    α^*  > max{1/e,p-2/p+2}

in what follows.



  
Stage 1
We will now determine k_1 such that α_k_1≥1/e.  We begin with a lemma.

 
For every α∈ (0,1), 

    H(α) > α^1-1/p.



We have

    H(α) 
       = p αμ(α)^p-1/ (p-1)μ(α)^p + α^p  
    = p αμ(α)^p-1/α-α^p/1-α + α^p  
    = p αμ(α)^p-1 (1-α) /α-α^p+1
       = p μ(α)^p-1 (1-α) / 1-α^p  
    = α^1-1/p g(α)^1-1/p/ h(α) ,

where

    g(α)    = 1-α^p-1/(p-1)(1-α) = 1/p-1∑_j=0^p-2α^j, 
    
    h(α)    = 1-α^p/p(1-α) = 1/p∑_j=0^p-1α^j.

Since 0 < h(α) < g(α) < 1 for every α∈ (0,1), it follows that

    g(α)^1-1/p/h(α) > g(α)/h(α) > 1.



Lemma <ref> implies

    α_k+1≥α_k^1-1/p

for every k, so 

    α_k ≥α^(1-1/p)^k.

Thus, we will have α_k_1≥ 1/e if α^(1-1/p)^k_1≥ 1/e, which means

    k_1 ≥loglog1/α/log (p/p-1).





  
Stage 2
As mentioned previously, k_2 is a constant independent of α and ϵ. 



  
Stage 3
We now determine k_3 such that 1-α_k_1+k_2+k_3/1+α_k_1+k_2+k_3≤ϵ.  By Lemma <ref>,

    1 - α_k+1/1+α_k+1
    = 1-H(α_k)/1+H(α_k)≤p/2( 1-α_k/1+α_k)^2

for k ≥ k_1+k_2.
In terms of δ_k := p/2( 1-α_k/1+α_k), we have δ_k+1≤δ_k^2, so

    δ_k_1+k_2+k≤δ_k_1+k_2^2^k≤( p/2( 1-α^*/1+α^*) )^2^k.

By (<ref>), p/2( 1-α^*/1+α^*) < 1, so we will have 1-α_k_1+k_2+k_3/1+α_k_1+k_2+k_3≤ϵ if 

    ( p/2( 1-α^*/1+α^*) )^2^k_3≤p/2ϵ,

i.e.

    k_3 ≥loglog2/ϵ p  - loglog2/p( 1+α^*/1-α^*) /log 2 .


Finally, by taking α=ϵ/2 we ensure that the error on [0,α^p] is bounded by ϵ (recall (<ref>)), so the error on [0,1] is bounded by ϵ. 

We illustrate the process in Figure <ref>, 
where we fix integers[p=31 is a somewhat arbitrary prime number, chosen in view of the number of days per month.] p and k, and numerically find the value of α∈ (0,1) and accordingly ϵ=1-α_k/1+α_k=2α such that with the (k,p,p-1)-composite rational approximant f_k the error is max_x∈[α^p,1]|x^1/p-f_k(x)|≤ϵ, achieved at x=1, and the error on [0,α^p] is bounded by ϵ. 
Observe that the maximum errors on [0,α^p] and [α^p,1] are 
 not equal but of the same order, suggesting the near optimality of our composite rational approximants.


Putting these inequalities together, we conclude that 

    k = loglog2/ϵ/log (p/p-1) + k_2 + loglog2/ϵ p/log 2

recursions are enough to yield accuracy ϵ, where k_2 is an integer satisfying

    k_2 - loglog2/p1+α^*/1-α^*/log 2≤k_2 ≤ k_2 - loglog2/p1+α^*/1-α^*/log 2 + 1.

Since k recursions translate into a rational function f_k of type (p^k-1,p^k-1-1), it follows that the degree n of the rational function f_k achieving accuracy ϵ is 

    n = p^loglog2/ϵ/log (p/p-1)+
    k_2+loglog2/ϵ p/log 2-1  .

We rewrite this to express the error with respect to the degree n. Taking the logarithm and absorbing the constant -1 into k_2, we get 

    log n = (loglog2/ϵ/log (p/p-1)+
    k_2+loglog2/ϵ p/log 2) log p ≤(loglog2/ϵ/log (p/p-1)+
    k_2+loglog2/ϵ/log 2) log p.

Hence,

    loglog2/ϵ≥log n-k_2log p/log p(1/log (p/p-1)+
    1/log 2) .

Thus, defining 

    c:= 1/log p(1/log (p/p-1)+
    1/log 2) = 
    log 2logp/p-1/log p log2p/p-1,

 we have 

    log2/ϵ≥(n/p^k_2)^c,

and therefore, writing b = 1/p^ck_2, we arrive at 

    ϵ≤ 2exp(-b n^c).

This bound holds when n is a sufficiently large power of p.  To handle the case in which n ∈ℕ is not a power of p, we note that ⌊log_p n ⌋ + 1 recursions yield a rational function of type (p^⌊log_p n ⌋, p^⌊log_p n ⌋-1), and for n large enough (n ≥ N, say) this function has error bounded above by

    2exp(-b (p^⌊log_p n ⌋)^c) ≤  2exp(-b p^-c n^c)  ≤exp(-(b p^-c-N^-clog 2) n^c).

Taking N large enough yields Theorem <ref> with b = b p^-c-N^-clog 2>0.

It is easy to see by comparing (<ref>) with (<ref>) that the same analysis, this time choosing α=ϵ rather than α=ϵ/2, also yields (<ref>) in Theorem <ref>. 


It remains to establish (<ref>). 
For this, we take α fixed and use a similar argument. In this case k_1,k_2 can both be regarded as constants independent of ϵ, since 

the error in the interval [0,α^p] is irrelevant. 
Therefore we write k:=k_1+k_2, and in place of (<ref>), the lowest degree n required for ϵ accuracy on S_p,α satisfies 
log n≤(k+loglog2/ϵ p/log 2) log p ≤(k+loglog1/ϵ/log 2) log p. Thus 
 defining 

    c:= log 2/log p (>c),

 we have 
log1/ϵ≥(n/p^k)^c,
and so setting b = 1/p^ck we obtain 
ϵ≤exp(-b n^c), as required. 
□

Note that for α≪ 1 we have k≈ k_1=
 loglog1/α/log (p/p-1), so b≈ (log1/α)^-C  for some C>0 depending only on p, so b scales like an inverse power of log1/α.



§ EXAMPLES
 
In Figure <ref> we illustrate our main result (<ref>) on approximation of x^1/p. For integers k=1,2,…,, we compute the error ϵ of the composite rational approximants as in Figure <ref>, 
and plot the errors against p^ck(≈ n^c) for p∈{2,5,31} in log-scale. The plots also show least-squares affine fits to the convergence data for each p. The fact that the affine fits closely trace the data suggests the exponent c in (<ref>) 
is sharp, especially for small values of p. For the p=31 plot, which ends early because computing further data was infeasible (note e.g.  that α^p< 10^-70 for k≥ 15), there is a slight bend in the convergence, which suggests that our c in (<ref>) might be a slight underestimate for large p.



Finally, Figure <ref> shows the error of the approximant g_k(z) to sect_p(z), which clearly exhibits equioscillation. Note how increasing k results in progressively smaller error (in log-scale), reflecting the double-exponential convergence. 
The error curves |g_k(z) - sect_p(z)| look identical on each of the segments [α,1]exp(2π ij/p) for j=0,…,p-1.







 §.§ Acknowledgment

We thank Alex Townsend, a discussion with whom inspired this work. 



abbrv


