The Fourth Operation

Use Ctrl + or – to enlarge or reduce text size.

The Fourth Operation
Introduction
Given whole numbers A and B we can add them together to get another whole number,

A + B

The definition of this operation is simply to count together what A and B count separately.

Then we can define multiplication as repeated addition:

A · B =	A + A + A + ...
	B times

And exponentiation as repeated multiplication:

A^{^B} =	A · A · A · ...
	B times

Those are the first, second, and third operations: addition, multiplication, and exponentiation. Presented in that way a prolonged series of operations comes to mind, each operation iterating the preceding one. The fourth operation would be (using a subscript to denote the second operand):

A_{_B} =	A^{^A}^{^{^{^A}}}^{^{^{^{^ททท}}}}
	B times

The series of operations can be defined recursively on the second operand.

Addition:

Given.

Multiplication:

A · 1 = A
A · (B + 1) = A + A · B

Exponentiation:

A^¹ = A
A^{^{B + 1}} = A · A^{^B}

The fourth operation:

A_₁ = A
A_{_{B + 1}} = A^{^A_{_B}}

The n^th operation, for n > 1:

A [n] 1 = A
A [n] (B + 1) = A [n – 1] ( A [n] B )

Historical note: A [n] B when A and B are whole numbers is the same as the function φ(A, B, n – 1) of Wilhelm Ackermann (1896 – 1962), used in a paper of 1928.

Thus we can define a series of operations between positive whole numbers. It is well known how to extend the definition of the first three operations to include all real numbers (except that in the case of exponentiation the base, A in A^{^B}, must be non-negative if the operation is to be well-defined and continuous as the exponent, B, varies).

The problem we investigate here is how to define the fourth operation, A_B or A raised to itself B times, for fractional B. What, for example, is A_½ or A raised to itself one half times? It might sound like a crazy question but no crazier than asking what is A multiplied by itself one half times. In that case the answer is the square root of A but the situation is not so simple with the fourth operation.

Later it will make some formulas simpler if we let A = e, the base of the natural logarithms.

e₀ = 1
e₁ = e = 2.71828…
e₂ = e^e = 15.15426...
e₃ = e^e^{^e} = 3814279.10476…

Again, what is

e_ฝ

e to itself one half times? And in general what is

e_s

for any real number s ?

... And what sort of nutball would ask such a question ?

Because it is there. Now Froggy, please stay out of the Algorithms section.

In order to solve the problem we will reformulate the fourth operation in terms of functional iteration, then we can use techniques of functional analysis to try to solve the problem. As well as our own ideas we will use some ideas from the paper by G. Szekeres “Fractional Iteration of Exponentially Growing Functions” (J. of the Australian Mathematical Soc., vol. 2, 1962), which relies on his “Regular iteration of real and complex functions” (Acta Math., vol. 100, 1958). Note that Szekeres’s point of view differs from ours in that his primary interest is (see below for definitions) in exp_s(x) as a function of x, where s is thought of as a parameter, whereas to us that function is only a stepping stone to exp_s(1) where s is the variable. Jumping ahead, note that Szekeres claims to have proved that the second coefficient in what we call Φ(x) is 1/2 but unless I’m mistaken he hasn’t done so or even made it plausible.

Sometimes we will write “exp” for the exponentiation of e function:

exp(x) = e^x

and place a subscript on exp to denote the level of its iteration. The functional iterates of exp(x) are closely related to e_s. First observe that

exp_₀(x) = x
exp_₁(x) = exp(x)
exp_₂(x) = exp(exp(x))
exp_₃(x) = exp(exp(exp(x)))

etc. Thus

exp _n (exp _m (x) ) = exp _n+m (x)

holds for all non-negative integers n and m. If it is to hold for all integers then exp_-1 must be the inverse of exp:

exp_{_-1}(x) = log(x)
exp_{_-2}(x) = log(log(x))

etc. where “log” is the natural logarithm function (in Calculus textbooks usually written ln ). Setting x = 1 gives us the fourth operation without the “socket”:

exp_n(1) = e^e^{^...}^{^{^e}}, n times = e_n

Thus our goal is to define the function

exp_s(x)

for all real numbers s, then

exp_s(1)

will answer our original question, what is e_s, that is, e raised to itself s times.

If that isn’t clear consider the third operation, exponentiation. It is iterated multiplication and can be obtained by iterating the function “mul”:

mul (x) = e x

Thus

mul (mul (x)) = e e x	or	mul₂(x) = e² x
mul (mul (mul (x))) = e e e x	or	mul₃(x) = e³ x

etc. In this case we know the continuous iterate of mul, it is

mul_s(x) = e^s x

so that

mul_s(1) = e^s

Iterates of the function of mul(x) where x is a sort of socket for the last multiplier gives a different way to look at exponentiation.

Groping for a Solution
First we will examine how the third operation, exponentiation, is defined when the exponent is a fraction and see if something analogous can be done for the fourth operation.

Unlike the operations of addition and multiplication, each of which has just one inverse (subtraction and division), exponentiation has two inverses because it is not commutative, the order of the operands matters. Given

A^{^B} = C

one operation is required to obtain the base A and another the exponent B. Taking a root gives the base

A = B^th root of C = ^{^B}√ C

and taking a logarithm gives the exponent:

B = logarithm (base A) of C = log_{_A}C

where contrary to usage later the subscript A does not denote iteration.

Consider the root inverse. Using facts about exponents with which I’ll assume you are familiar, given that we have defined exponentiation when the exponent is an integer we can use the root operation to define it when the exponent is any rational number. In what follows let m and n be integers, n positive. Then let

C^1/n = ⁿ√ C

and

C^m/n = ⁿ√ C ^m

This comports with the definition of exponentiation when m/n is an integer and otherwise obeys the expected behavior of exponentiation. (If you investigate you will see it does this because multiplication is commutative and associative, multiplication doesn’t depend on the order or grouping of its operands.) Then by insisting that exponentiation be continuous we can define it for real s by taking any sequence of rational numbers r_i that converges to s as i → ∞ (here the subscript on r is an index, not the fourth operation):

lim r_i = s
i → ∞

and defining

C^{^s} = lim C ^{^r _i}
i → ∞

(one would have to prove that this is well defined, that it does not depend on the sequence).

The trouble with generalizing this procedure (using root, the inverse for the first operand) is that, to repeat, it uses the fact that multiplication is commutative and associative. These two properties give a certain homogeneity to multiplication which allows us to turn a fraction 1/n into an integral root, and a fraction m/n into an integral power of that root. This will not work for the fourth operation because exponentiation is so different from multiplication.

Now consider the other inverse of exponentiation, the logarithm. Though its basic definition is that it is the exponent there is an interesting theorem which allows us to compute it more directly. Let

f (x) = ∫_₁^{^x} 1/t dt

Then for all real x and y

f (x y) = f (x) + f (y)

(On the right side, in the second integral make the change of variables u = xt .) Thus for any whole number p

f (x^p) = p f (x)

Also

q f (x^1/q) = f (x^1/q^{^q}) = f (x^q^{^1/q}) = f (x)

so that

f (x^1/q) = 1/q f (x)

Thus for any rational number r

f (x^r ) = r f (x)

and since f is continuous (indeed differentiable) this holds for any real number r. Finally, let e be the solution to f (x) = 1 (there is one and only one solution because f (1) = 0 and f (x) increases to infinity). Then

f (e^r ) = r

Thus f (x) is the logarithm for exponentiation using the base e and we have a very interesting formula for it:

log x = ∫_₁^{^x} 1/t dt

And we have e^x as well because it is the inverse of log x.

Note that the derivative of log x is:

log^’x = 1/x

a simple function independent of log x and of e^x. Thus we can use 1/x to determine log x. This contrasts with the derivative of e^x , which is e^x again.

Though the second procedure (using logarithm, the inverse for the second operand) doesn’t generalize to higher operations, for the same reason the first doesn’t, it does furnish an important clue: if one seeks to continuously iterate a function – any function, not necessarily exp – then instead of trying to find it directly, look for the derivative of its inverse – the inverse that gives the “higher” operand. In other words, try to define the derivative of the generalized “logarithm” that counts the number of times the function was applied in its argument.

Uniqueness
The hard part of the problem of defining the fourth operation e_s isn’t finding a definition it is finding a natural, best definition. There are any number of “low quality” solutions: in the interval 0 ≤ s < 1 set

e_s = any function of s that is one when s is zero

then use the fact that

e_s+1 = e^e_s

to translate the values at 0 ≤ s < 1 to any higher unit interval. In effect, for s ≥ 1 find its integral part n and fractional part t and set

e_s = exp_n(e_t)

What is wanted is a condition on our definition of e_s that makes it the unique best definition.

The problem is analogous to defining the factorial function

n ! = n (n – 1) (n – 2) ... 1

for fractional n. One could define 0 ! = 1 and s ! as anything when s is between 0 and 1, then use

(s + 1) ! = s ! (s + 1)

to define it for other values of s. However one particular definition of s! (we won’t give the details) has two recommendations: (1) it arises by a simple and natural argument, (2) it is the only function defined for s ≥ 0 such that 0 ! = 1, s ! = s (s – 1) !, and that is log convex. It is also the only such function that is asymptotically equal to a logarithmico-exponential function. As Szekeres points out in the first reference above, this analogy suggests that in the search for a best solution of fractional functional iterates it is “better to concentrate on the real variable and asymptotic properties than on the complex-analytic character of the solution.” We seek a real variable condition that will make the fourth operation the fourth operation.

A Related Problem
Finding a continuous iterate of e^x that is natural and arguably the best is a difficult problem to attack directly, so instead at first we will try to find such a continuous iterate for a function not much different

e(x) = e^x – 1

This function has a fixed point, and only one fixed point, at zero:

e(0) = 0

It is much easier to find a continuous iterate of a function with a unique fixed point.

... I knew a man who dropped his latchkey at the door of his house and began looking for it underneath the street light. I asked him why he was looking for it there when he knew very well it was over by the door. He replied “Because there isn’t any light over there.” LOL.

In our case, Froggy, we might find a key to the key where the light is. Perhaps if we solve the continuous iteration problem for e(x) we will be able to use it to solve the problem for exp(x).

Our choice of notation can be confusing in that sometimes e stands for the number 2.718 ... and sometimes the function e(x). e_s(x) is a function and the function it iterates is exp(x) – 1, yet if we were to write the function without its argument, e_s, it would look like the fourth operation — e^(e^(...^e)...) s times — rather than the s^th iterate of the function e(x) = exp(x) – 1. To avoid confusion we will never do that, that is, we will always include at least a placeholder dummy variable when referring to the function e(x) and its iterates.

Finding the “logarithm” for e_s(1) is simpler than finding e_s(1). We mean logarithm in a general sense, as counting the number of times something has been done. The ordinary logarithm of Calculus counts the number of factors of e in a product, for example

log(e e e e) = 4

and in general

log(e^s) = s

where s can be fractional as well as whole. The crux of the definition, see the discussion of mul(x) = e x above, is that

log(mul (x) ) = log(x) + 1

Facts About Functions
When speaking of a general function we shall always assume it is strictly increasing and smooth (infinitely differentiable). In general the “logarithm” for iterates of a function f (x) is a function A(x) defined for x > 0 such that

A( f (x) ) = A(x) + 1

and

A is strictly increasing and smooth
A(1) = 0

The functional equation is called Abel’s equation, after Niels Henrik Abel (1802 – 1829). A(x) counts the number of times f (x) is used in x, it is a sort of “functional logarithm.” (Some authors don’t include A(1) = 0 in the definition, and when it is posited call the resulting A(x) a “normalized” Abel function.)

You can see from Abel’s equation that:

A(f₂(x)) = A( f ( f (x) ) ) = A( f (x) ) + 1 = A(x) + 2

and in general

A( f_n(x) ) = A(x) + n

Thus if we had such a function A we could define f_s for real s as the solution to

A( f_s(x) ) = A(x) + s

that is,

f_s (x) = A _-1( A(x) + s )

where A _-1 is the inverse of A. Such an f_s (x) would have the properties

f_₀ (x) = x
f_s (f_r (x)) = f_s + r (x)

Note that

f_s (1) = A _-1(s)

That works in reverse. If f_s(x) is a continuum of iterates for f (x), that is, f_₀(x) = x and f_s (f_r (x)) = f_s + r (x) , then f_s can be obtained from the inverse of the function

B(s) = f_s(1)

as above. First replace s with s + r :

B(s + r) = f_s + r (1) = f_s(f_r(1))

Given x let r = B _-1(x). Then

B(B _-1(x) + s) = f_s(x)

Let A = B _-1 and we have our original formula.

If we knew nothing more about A than Abel’s equation then A would not be unique. If (given the function f ) A is one solution of Abel’s equation then so would be B(x) = A(x) + g(A(x)) for any periodic function g with period 1 and g(0) = 0. This is easily seen by substituting B into Abel’s equation. (Abel showed that this describes all solutions of Abel’s equation.)

As might be expected from the formula for the logarithm for exponentiation the derivative of A(x) might be useful. Call

a(x) = A^’(x)

a “derived Abel function” of f. Let

Φ(x) = 1 / a(x)

Then because of the condition A(1) = 0 we have

A (x) = ∫_₁^{^x} 1 / Φ(t) dt

At this point the equation looks contrived but if we can determine Φ(x) we solve the iteration problem for f(x). It turns out we will be able to do it when f (x) = e^x – 1 .

If we take the derivative of each side of Abel’s equation, after some manipulation we get

Φ(e(x)) = f ^’(x) Φ(x)

Call this Schr๖der’s equation for f (x) , after Ernst Schr๖der (1841 – 1902).

Although we won’t need it here, it’s interesting that if we let

f_s(x) = Φ _-1( f ^’(s x) Φ(x) )

then if f ^’(0) = 1 the above is a system of continuous iterates of f (x).

Recall that if a function f (x) is smooth there is a power series that might not converge but at least will be asymptotic to the function at 0 :

f (x) = f (0) + f ^’(0) x + 1/2 f ^’’(0) x² + 1/6 f ^’’’(x) x³ + ...
as x → 0

We will explain what asymptotic means in the next section. The general term is

1 / i ! f ^⁽ⁱ⁾(0) x^ⁱ

where f ^⁽ⁱ⁾(x) denotes the i^th derivative of f (x).

If the function is not only smooth but real analytic the series converges and we needn’t specify “as x goes to zero.” Some authors write ≈ instead of = when the series diverges or isn’t known to converge. We shall let “as x goes to zero” serve to indicate that the series is asymptotic but might not converge.

Asymptotic Power Series
Suppose we are given, out of the blue, an infinite series

a₀ + a₁(x – x₀) + a₂(x – x₀)² + a₃(x – x₀)³ + ...

Regard this as a “formal” series, effectively just a sequence of coefficients provenance unknown.

Given a particular x we could try to sum the series. The series might converge, that is, if we compute the partial sums

S_n(x) = a₀ + a₁(x – x₀) + a₂(x – x₀)² + ... a_n(x – x₀)ⁿ

then

lim S_n(x) = [some number] as n → ∞

for the particular x. In such cases we can think of the series as a function defined where it converges. The series would stand alone and could be used to define the function it converges to.

The series always converges at x = x₀, but except for that value it might not converge. If it doesn’t converge anywhere besides at x₀ it is useless, by itself.

Scrap the above series and suppose we are given a function F(x) instead of a series. We might not know how to compute F(x) but suppose we have an abstract definition, know what its properties are. And suppose we know it is smooth at x₀, any number of its derivatives exist there. Then we can form the Taylor expansion of F at x₀

F(x₀) + F^’(x₀) (x – x₀) + 1/2 F^’’(x₀) (x – x₀)² + ... 1 / i ! F⁽ⁱ⁾(x₀)(x – x₀)ⁱ + ...

As before we regard this as a formal series that is not necessarily summable, in effect just a sequence of coefficients. For a finite number of terms we can compute the partial sum

S_n(x) = F(x₀) + F^’(x₀) (x – x₀) + 1/2 F^’’(x₀) (x – x₀)² + ... 1 / i ! F⁽ⁿ⁾(x₀)(x – x₀)ⁿ

The sequence might diverge, that is, if we fix x ≠ x₀ and let n → ∞ the limit of S_n(x) might not exist. But even so the series tells something about F near x₀. Instead of fixing x and letting n get arbitrarily large, if we fix n and let x get arbitrarily close to x₀, then by definition of the derivative S_n(x) will approximate F(x).

But what good is that? We already know the series gives the correct value when x = x₀, what use is an approximation for x near x₀ if what we want to know is F(x) for x well away from x₀? No use at all, unless – and this is the point – the function possesses a property that allows us to put its value at x in terms of its value at points near x₀.

If the property allows us to approach x₀ arbitrarily close then the asymptotic series, despite the fact that it doesn’t converge, can be used to compute F(x) to any desired accuracy. We will see an example of this in the next section.

To repeat, with a convergent series, given any x (in the domain of convergence) the series converges as n → ∞. With an asymptotic series, given any n the truncated series converges as x → x₀. A convergent series defines a function, an asymptotic series by itself does not – however it can be used to compute a function whose value near x₀ determine its value far from x₀.

Again, convergence is concerned with S_n(x) for fixed x as n → ∞ whereas being asymptotic is concerned with S_n(x) for fixed n as x → x₀ when one knows the series is the Taylor expansion of a function that we know other properties of besides the Taylor coefficients. Convergence depends on only the coefficients, asymptotic depends on the function as well as the coefficients.

We have been using “asymptotic” as if the Taylor series diverges but since we might not know that in advance it is useful to subsume convergent series under asymptotic series, considering that the definition of convergence generally includes the condition for asymptotic.

To Proceed
As before, let

e(x) = e^x – 1

It has a fixed point at 0. Its power series expansion about 0 is

e (x) = x + 1/2 x² + 1/6 x³ + ...

Let A(x) be an Abel function for e(x),

A(e(x)) = A(x) + 1

and a(x) the corresponding derived Abel function for e(x),

a(x) = A^’(x)

Let

Φ(x) = 1 / a(x)

Schr๖der’s equation for e(x) is

Φ(e(x)) = e^’(x) Φ(x)

and using the fact that e^’(x) = e^x

Φ(e(x)) = e^x Φ(x)

We will try to find the Taylor series for Φ(x) at 0. We can expect only that it converges asymptotically, nonetheless, as we shall see later, knowing Φ(x) when x is near 0 will allow us to compute it for any x .

Φ(x) = c_₀ + c_₁ x + c_₂ x² + c_₃ x³ + ... as x → 0⁺

It is fairly easy to show that both c_₀ and c_₁ must be 0. Take the derivative of both sides of Schr๖der’s equation to eventually obtain

Φ^’(e(x)) = Φ^’(x) + Φ(x)

Then by considering the particular case x = 0 we can conclude that

Φ(0) = 0

Then, returning to the last equation in x , take the derivative of each side again, again evaluate at x = 0 and we can conclude that

Φ^’(0) = 0

Continuing in the same vein yields nothing about Φ^’’(0) but at least now we know that the series begins

Φ(x) = c_₂ x² + c_₃ x³ + c_₄ x⁴ + ... as x → 0⁺

where we have yet to determine the remaining coefficients c_₂, c_₃, c_₄, ...

To determine them we will plug the known power series for e(x) and the unknown one for Φ(x) into Schr๖der’s equation, then equate coefficients of the same power of x and solve for c_₂, c_₃, etc. one after the other. The trouble is that when we do it there is one degree of freedom. That is, we need to know at least one of the c_{_i} in order to solve for the others. If we knew c_₂ we could solve for c_₃, then for c_₄, etc. We will explain in detail when we have figured out what c_₂ should be and actually perform the procedure.

Finding c_₂ in the Taylor series for Φ(x) will take a bit of work. What we will do is find the beginning of the Taylor series for e_s(x) in two different ways, first thinking of x as the variable then s. In the first case the coefficients will be functions of s, in the second functions of x. Looking ahead, in the first case we will be able to determine those functions, in the second they will involve Φ(x). By equating the series and comparing the same powers of s, we can make a plausible choice for c_₂ in Φ(x). However, in the last analysis it looks like it will have to be a condition, an imposed value, rather than something we can deduce. But given this one value everything else is determined. Φ(x) will be unique given c_₂ .

To begin, let’s see how far we can go finding a power series for e_s(x). We can try to find one in powers of x and in powers of s.

First consider e_s(x) as a function of x with s as a parameter. Consider its Taylor expansion at x = 0

e_s(x) = c_₀ + c_₁ x + c_₂ x² + c_₃ x³ + ... as x → 0

The coefficients c_{_i} are functions of s , namely the i^th derivative of e_s(x) with respect to x divided by i ! evaluated at x = 0 .

Since e₀(x) = x for all x, including when x = 0, it must be that c_₀ = 0 and we can write

e_s(x) = c_₁ x + c_₂ x² + c_₃ x³ + ...

Now c_₁ is the first derivative of e_s(x) evaluated at x = 0. To try finding the first derivative we use a trick: before taking the derivative break e_s(x) into e(e_s–1(x) ) then use the chain rule to find a recursive relation for the derivative. In what follows the apostrophe indicates differentiation with respect to x.

e_s(x) = e(e_s–1(x))
e^’_s (x) = e^’(e_s–1(x)) e^’_s–1 (x)

Then since e_s–1(0) = 0 and e^’ (0) = 1

e^’_s (0) = e^’_s–1 (0)

Repeating that formula within itself on the right we get down to

e^’_s (0) = e^’_t (0)

where 0 ≤ t < s is the fractional part of s, that is, the number left over after subtracting the largest whole number less than s from s.

We have failed to determine e^’_s (0) in general. However in the particular cases when s is a whole number, that is when t = 0, we have determined it:

e^’_n (0) = 1

for every whole number n. We shall impose the natural, simplest, condition on e_s(x) that generalizes that formula:

Condition #1
e^’_s (0) = 1

With that condition c_₁ = 1 independent of s. Thus the power series for e_s(x) begins:

e_s(x) = x + c_₂ x² + c_₃ x³ + ...

Now try the same trick with the second derivative. If we can find it c_₂ will be it divided by 2. For typographical reasons in what follows let D as well as apostrophe denote d/dx, differentiation by x, and let D^² as well as double apostrophe denote differentiation twice.

D^²[ e_s(x) ] = D[ e^’(e_s–1 (x) ) e^’_s–1 (x) ]
= e^’’(e_s–1 (x) ) e^’_s–1 (x) e^’_s–1 (x) + e^’_s–1 (x) e^’’_s–1 (x)

Setting x = 0 and using the fact that e_{s – 1}(0) = 0 and e^’_{s – 1} (0) = 1 and e^’’ (0) = 1, then simplifying, we get

e^’’_s (0) = 1 + e^’’_s–1 (0)

Repeating that formula within itself on the right

e^’’_s (0) = 1 + 1 + e^’’_s–2 (0)

etc. so that

e^’’_s (0) = n + e^’’_t (0)

where n + t = s, n being the largest whole number ≤ s and 0 ≤ t < s .

We have failed to determine e^’’_s (0) in general. However in the particular cases when s is an integer, that is, t = 0, n = s, we have determined it:

e^’’_n (0) = n

We impose the natural, simplest, condition on e_s(x) that generalizes that formula:

Condition #2
e^’’_s (0) = s

Thus we have the power series for e_s(x) beginning:

e_s(x) = x + 1/2 s x² + c_₃ x³ + ... as x → 0⁺

We will not need any more conditions. The rest of the coefficients can be determined by breaking e_s+1(x) two ways

e(e_s(x)) = e_s(e(x))

and plugging in the known power series for e(x) and what we know of the power series for e_s(x) and equating coefficients of like powers of x. It helps to have a symbolic algebra calculator such as PARI/GP by the PARI Group of the Institute of Mathematics at the University of Bordeaux. The next three coefficients are

c_₃ = 1/4 s² – 1/12 s
c_₄ = 1/8 s³ – 5/48 s² + 1/48 s
c_₅ = 1/16 s⁴ – 13/144 s³ + 1/24 s² – 1/180 s

so that

e_s(x)^²

^² =

x + 1/2 s x² + (1/4 s² – 1/12 s) x³ + (1/8 s³ – 5/48 s² + 1/48 s) x⁴ + (1/16 s⁴ – 13/144 s³ + 1/24 s² – 1/180 s) x⁵ + ...

as x → 0

For any finite number of terms we can rearrange in powers of s.

e_s(x)^²

^² =

x + (1/2 x² – 1/12 x³ + 1/48 x⁴ – 1/180 x⁵ + ...) s + (1/4 x³ – 5/48 x⁴ + 1/24 x⁵ + ...) s² + (1/8 x⁴ – 13/144 x⁵ + ...) s³ + (1/16 x⁵ + ...) s⁴ + ...

still as x goes to 0, not s .

Now consider e_s(x) as a function of s with x as a parameter. Let its Taylor expansion at s = 0 be

e_s(x) = c_₀ + c_₁ s + c_₂ s² + c_₃ s³ + ... as s → 0

The coefficients c_{_i} are functions of x , namely the i^th derivative of e_s(x) with respect to s evaluated at s = 0 divided by i ! .

Since e₀(x) = x right away we know c_₀ = x .

Now c_₁ is the derivative of e_s(x) with respect to s evaluated at s = 0. To find that consider an Abel function for e(x), and we say “an” rather than “the” because at this point we do not know if it is unique

e_s(x) = A_{_-1}( A(x) + s )

where A_{_-1} is the inverse of A. We will differentiate each side with respect to s, using the fact that the derivative of the inverse of a function is the reciprocal of its derivative evaluated at the inverse. We will also use the fact that the derivative of A(x) + s with respect to s is 1 . In the next few equations an apostrophe indicates differentiation with respect to s.

e_s^’ (x) = 1 / A^’(A_{_-1}(A(x) + s) )
= 1 / A^’(e_s(x))
= 1 / a(e_s(x) )

Thus since e₀(x) = x

e₀^’ (x) = 1 / a(x) = Φ(x)

So the power series in s begins

e_s(x) = x + Φ(x) s + c_₂ s² + c_₃ s³ + ... as s → 0

It’s easy to get all the higher derivatives with respect to s from the formula for the first derivative, e_s^’ (x) . You find that they form a sequence of odd powers of Φ(x) times a constant (depending on i ) that alternates in sign. By way of a footnote (we won’t need it) the exact formula for the i^th derivative of e_s(x) with respect to s when i > 1 is

– (–1)^ⁱ (2i–3)!! Φ^{^2i–1}(e_s(x))

where !! indicates the double factorial: N!! = N (N – 2) (N – 4) (N – 6) ... and the last factor is 1 or 2 depending on whether N is odd or even; in the case here it’s always odd. That i^th derivative divided by i! and evaluated at s = 0 gives c_{_i} for i > 1

c_{_i} = – (–1)^ⁱ (2i–3)!! / i! Φ^{^2i–1}(x)

We are trying to find Φ(x) . We proved that its Taylor expansion begins

Φ(x) = c_₂ x² + c_₃ x³ + c_₄ x⁴ + ... as x → 0⁺

and if we can determine just c_₂ we can determine the others. Now we found formulas for e_s(x) in two ways and the second involved Φ(x).

e_s(x) = x + 1/2 s x² + [higher powers of s and x] as x → 0
e_s(x) = x + Φ(x) s + ... as s → 0

It would solve our problem if we could equate coefficients of the first power of s and conclude that

Φ(x) = 1/2 x² – 1/12 x³ + 1/48 x⁴ + ... as x → 0⁺

but it isn’t obvious that they must be equal. For one thing the first series is as x goes to 0 and the second as s goes to 0. Still it makes plausible that c₂ in the power expansion of Φ(x) is 1/2. Furthermore, equating coefficients of the first power of x rather than of s results in the very same value. Therefore we impose the following condition, which supersedes our previous conditions #1 and #2 because alone it will allow us to determine Φ(x).

Condition
The second order coefficient in the power expansion of Φ(x) is 1/2.

Thus we have

Φ(x) = 1/2 x² + c_₃ x³ + c_₄ x⁴ + ... as x → 0⁺

Now to find the remaining coefficients. Consider Schr๖der’s equation

Φ(e(x)) = e^x Φ(x)

We know the power series for e^x

e^x = 1 + x + 1/2 x² + 1/6 x³ + ...

and for e(x)

e (x) = x + 1/2 x² + 1/6 x³ + ...

Begin by thinking of Φ(x) as

Φ(x) = 1/2 x² + c x³ + [higher powers of x]

Then plug these power series into Schr๖der’s equation. (You might want to use a symbolic algebra computer program such as PARI/GP mentioned earlier.) The left side is

1/2 x² + (c + 1/2) x³ + (3/2 c + 7/24) x⁴ + [higher powers of x]

and the right

1/2 x² + (c + 1/2) x³ + (c + 1/4) x⁴ + [higher powers of x]

Thus, focusing on the fourth order term,

3/2 c + 7/24 = c + 1/4

we have

c = – 1/12

Then think of Φ(x) as

Φ(x) = 1/2 x² – 1/12 x³ + c x⁴ + [higher powers of x]

and again plug into Schr๖der’s equation. As before the first few coefficients are tautological but equating the ones for x⁵

2 c + 1/48 = c + 1/24

and solving for c

c = 1/48

Continuing in this manner we can determine the coefficients one after the other.

Φ(x)²

1/2 x² – 1/12 x³ + 1/48 x⁴ – 1/180 x⁵ + 11/8640 x⁶ – 1/6720 x⁷ – 11/241920 x⁸ ...

as x → 0⁺

Remarkably, this appears to be the first degree coefficient of the power series of e_s(x) in s. Positing that the second order term in Φ(x) be the same as the corresponding term in e_s(x) seems to entail that all of them are the same. I say “seems” because we can only expand to a finite number of terms and so far the conjecture holds. If we could prove the conjecture then we would not need our condition on “c₂.” In any case, expanding Φ(x) directly is much easier than expanding e_s(x); it has the advantage of being a function of only one variable.

The coefficients continue (order 9 to 13): 29/14551520, 493/43545600, –2711/239500800, –6203/3592512000, 2636317/373621248000.

Now that we have an asymptotic power series for Φ(x) as x → 0⁺ we can use it, together with known properties of the function Φ(x), to calculate Φ(x) for any x .

First note that the inverse e_{_–1}(x) of e(x) = exp(x) – 1 is

e_{_–1}(x) = log(1 + x)

Repeated application of this function takes a number arbitrarily close to zero:

lim e_{_– n}(x) = 0
n → ∞

Now suppose we are given x = x_₀ > 0. Here a subscript on x is a sequential index, not the fourth operation. Then the sequence

x_₀
x_₁ = e_{_–1}(x_₀)
x_₂ = e_{_–2}(x_₀)
...
x_{_n} = e_{_– n}(x_₀)

can be continued far enough (n can be made large enough) so that x_{_n} is as small as desired.

Now apply Schr๖der’s equation

Φ(e(x)) = e^x Φ(x)

to x_₁

Φ(x_₀) = Φ(e(x_₁)) = e^x_₁ Φ(x_₁)

then to x_₂

Φ(x_₁) = Φ(e(x_₂)) = e^x_₂ Φ(x_₂)

Thus

Φ(x_₀) = e^x_₁ e^x_₂ Φ(x_₂)

And so on through x_{_n}

Φ(x_₀) = e^x_₁ e^x_₂ ... e^x_{_n} Φ(x_{_n})

Since x_{_i+1} = e_{_–1}(x_{_i}) = log(1 + x_{_i}) we have

Φ(x_₀) = (1 + x_₀) (1 + x_₁) ... (1 + x_{_n–1}) Φ(x_{_n})

We can use this equation to, in a manner of speaking, drop any x = x_₀ > 0 down to near zero. Thus with our asymptotic expansion of Φ(x) at 0 we can compute Φ(x) for any x.

Having Φ(x) we can compute the “functional logarithm” for e(x) from

A (x) = ∫_₁^{^x} 1 / Φ(t) dt

Then the functional iterates of e(x) are

e_s(x) = A _{_– 1}(A(x) + s)

Note that

e_s(1) = A _{_– 1}(s)

And because e_s(1) = e_n(e_{s – n}(1) ) we have

A _{_– 1}(s) = e_n(A _{_– 1}(s – n) )

so if we know A _{_– 1}(s) for 0 ≤ s < 1 we know it for all s.
Φ(x) is unique (given the “c₂ condition”) because we have calculated what it must be. Thus also A(x) and its inverse A _{_– 1}(x) are unique.

Uniqueness Again
Since the uniqueness of Φ(x) is so important we shall prove it again without resorting to an asymptotic expansion, simplifying an argument of Szekeres.

Let Φ(x) be the reciprocal of a derived Abel function for e(x). We shall prove that if Φ(0) = Φ^’(0) = 0, and Φ^’’(0) exists and ≠ 0 then Φ(x) is unique.

The condition says

lim Φ(x) / x^² exists and ≠ 0
x → 0

Call this constant c .

Suppose there was another derived Abel function for e(x) whose reciprocal, call it Φ^*(x), satisfied our condition. We want to show it must be that Φ^* = Φ.

Choose one x = x_₀ and consider the sequence

y_{_n} = e_{_– n}(x_₀)

It goes to 0 (see earlier text). We will use Schr๖der’s equation for each Φ to express y_{_n} in terms of x_₀. It’s easy to see from the equation, by substituting e_{_– n}(x) for x, that

Φ(e_{_– 1}(x)) = e^–x Φ(x)

Thus

Φ(e_{_– n}(x)) = e^–n Φ(x)

So for each Φ

Φ(y_{_n}) = e^–n Φ(x_₀)
Φ*(y_{_n}) = e^–n Φ*(x_₀)

Thus their ratio doesn’t depend on n. Call the ratio r.

r = Φ*(y_{_n}) / Φ(y_{_n}) = Φ*(x_₀) / Φ(x_₀)

so that

Φ*(y_{_n}) = r Φ(y_{_n})

At this point it looks like r depends on x_₀.

By the condition we know

lim Φ(y_{_n}) / y_{_n}² = c
n → ∞

where c ≠ 0. Hence

lim Φ*(y_{_n}) / y_{_n}² = r / c
n → ∞

The left side is a constant so r is the product of two constants. Thus r is constant and independent of x_₀, and

Φ*(x) / Φ(x) = r

Φ*(x) = r Φ(x)

for every x.

Now recall that for each Φ and all x

∫_₁^{^x} 1 / Φ(t) dt = A (x)

Letting x = e(1)

∫_₁^{^e(1)} 1 / Φ(t) dt = A (e(1)) = 1

Thus

∫_₁^{^e(1)} 1 / Φ*(t) dt = 1 / r ∫_₁^{^e(1)} 1 / Φ(t) dt

1 = 1 / r

Thus r = 1 and Φ*(x) = Φ(x) for all x.

The Original Problem
We began by asking what is the number e to itself a fractional number of times. We reformulated the problem as finding the fractional iterates of the function exp(x). Then instead of solving that problem we found the fractional iterates of a somewhat different function e(x) = exp(x) – 1.

Knowing how to compute the fractional iterates of e(x) will enable us to compute the fractional iterates of exp(x). The idea is to use integral iterations of exp(x), which we know already, then introduce a fractional iteration of e(x) which will differ a little from the desired fractional iteration of exp, then take an iterated log so the error is reduced and the original integral iteration of exp is undone. The error goes to zero as the number of integral iterations increases. We now give the details.

We need define exp_s(x) only for 0 ≤ s < 1 because exp_s(x) for larger values of s can then be found by repeated application of

exp_s+1(x) = exp(exp_s(x) )

and for smaller values by repeated application of

exp_s–1(x) = exp_–1(exp_s(x) ) = log(exp_s(x) )

So in groping for a definition we will at first suppose 0 ≤ s < 1. At this point all we impose on exp_s(x) is that it be an increasing function of s and of x, and

exp_r(exp_s(x) ) = exp_r+s(x)

for r and s fractional just as when they are integral.

Recall that e(x) = exp(x) – 1. Consider the identities

exp₀(x) = e₀(exp₀(x) ) + 0

and

exp₁(x) = e₁(exp₀(x) ) + 1

The first just says x = x and the second exp(x) = e(x) + 1. Since exp_s(x) and e_s(x) are both increasing functions of s the following condition on exp_s(x) is reasonable.

Assumption
If 0 ≤ s < 1 then
exp_s(x) = e_s(x) + ε
where 0 ≤ ε < 1.

Now let n be a whole number and consider

exp_s(x) = exp_–n(exp_n+s(x) ) = exp_–n(exp_s(exp_n(x) ) )

Since exp_–n is log iterated n times,

exp_s(x) = log_n(exp_s(exp_n(x) ) )

The idea is to replace the exp_s() on the right with e_s() and let the iterated exp / log wash out the difference as n grows large.

By our assumption above, replacing x with exp_n(x)

exp_s(x) = log_n(e_s(exp_n(x) ) + ε)

where 0 ≤ ε < 1 depends on n and x as well as s. We want to prove that the difference between the right side and

log_n(e_s(exp_n(x) ) )

goes to 0 as n → ∞.

This is a consequence of the increasingly slow growth of log_n(x) as either n or x increases. First note that

log^’x = 1 / x

for x > 0, and in general

log_n^’x = 1 / x log x log_₂x log_₃x ... log_n–1x

for x large enough that the arguments of all the logarithms exceed zero. Define the sequence of numbers

a_n = exp_n –1(x)

Then

log_n^’a_n = 1 / exp_n–1x exp_n–2x exp_n–3x ... exp₀x → 1

as n → ∞. And since log_n^’x is a decreasing function of x the above limit is true for any sequence whose terms are greater than a_n, such as e_s(exp_n(x)) > e_s(exp_n-1(x)) ≥ exp_n-1(x).

Letting y = e_s(exp_n(x)) we have

0 < log_n(y + ε) – log_n(y) < ε log_n^’(y)

and the right side goes to zero as n goes to infinity.

Thus

exp_s(x) = lim log_n(e_s(exp_n(x) ) )
n → ∞

In practice the limit converges very rapidly. So rapidly in fact that n = 2 or 3 is enough for eight significant figures.

Here is how much exp_s(x) differs from e_s(x) when s is between 0 and 1. To three decimal places:

x =	1	2	3	4	5	6	7	8	9	10
s	ε	ε	ε	ε	ε	ε	ε	ε	ε	ε
0.00	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
0.25	0.173	0.114	0.064	0.031	0.013	0.006	0.002	0.001	0.000⁺	0.000⁺
0.50	0.374	0.261	0.150	0.075	0.036	0.017	0.008	0.003	0.002	0.001
0.75	0.633	0.498	0.340	0.217	0.134	0.080	0.047	0.027	0.015	0.008
1.00	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000

Naturally if s = 1 the difference is always one but if s < 1 not only is the difference less than one it looks like it goes to zero, and rapidly, with increasing x .

exp_s(x) thus defined is a continuum of iterates for exp(x). Even for finite n we have

log_n⚬ e₀ ⚬ exp_n = log_n ⚬ I ⚬ exp_n = log_n ⚬ exp_n = I

and in

log_m⚬ e_r ⚬ exp_m ⚬ log_n⚬ e_s ⚬ exp_n

since the limit on m and on n exists and is unique we can choose n = m so in the limit on m the above is

log_m⚬ e_r ⚬ e_s ⚬ exp_m = log_m⚬ e_r+s ⚬ exp_m

Now to answer the question we started with, what is e to itself one half times, we compute exp_ฝ(1) using the above limit formula.

n	formula
0	1.2710274
1	1.6100969
2	1.6451508
3	1.6451508

the last two accurate to seven decimal places. As a check that this great web of arithmetic hangs together we compute exp_ฝ(exp_ฝ(1) ). By n = 2 it has stabilized at seven decimal places to 2.7182818, which is as it should be exp(1) = e, the base we are using.

The following table shows how fast the fourth operation grows as the second operand increases. The numbers in the right column are to three decimal places. The zeros at the end of the higher numbers are just placeholders due to the limited precision of the computer.

s e_s
0.000     1.000
0.125     1.141
0.250     1.293
0.375     1.460
0.500     1.645
0.625     1.855
0.750     2.097
0.875     2.380
1.000     2.718
1.125     3.130
1.250     3.644
1.375     4.305
1.500     5.182
1.625     6.392
1.750     8.141
1.875     10.807
2.000     15.154
2.125     22.879
2.250     38.261
2.375     74.059
2.500     178.001
2.625     597.197
2.750     3431.358
2.875     49362.812
3.000     3814279.105
3.125     8632298240.322
3.250     41348087974198336.000
3.375     145673789637930926000000000000000.000
3.500     201844519565335862000000000000000000000000000000000000000000000000000000000000.000

Here is a graph of y = e_x with the scale of the y axis shrunk in order to show more of the graph.

Features to note:

e_x is undefined for x ≤ –2. (It would involve the logarithm of a non-positive number.)

As x → –2⁺, e_x → –∞.

When x = –1, e_x = 0.

When x = 0, e_x = 1.

When x = 1, e_x = 2.71828 ... = e.

There is a point of inflexion between –1 and 0. Numerical analysis (specifically, examining the second differences) shows it is at about x = –.5838, y = .1710 .

By definition a graph is fairly straight near a point of inflexion. This graph is fairly straight throughout the entire unit interval from x = –1 to 0. Throughout that range e_x differs from x + 1 by less than .01 .

Is it Unique?
The question of uniqueness bedevils attempts to define the fourth operation. Earlier we showed that, given some justified assumptions, e_s(x) is the unique continuum of iterates of e(x). Then we defined exp_s(x) in terms of e_s(x). It isn’t clear though that this exp_s(x) is unique. Without proof, there might be another way to define it, using e(x) or not.

Other Terminology and Notation
The fourth operation is known by various names: tetration, super-exponentiation, hyperpowers.

Instead of our subscript notation e_s some authors use the pre-superscript notation ^se (Hans Maurer [1868 – 1945] used it in a paper of 1901). Other notation is Donald Knuth’s up arrow notation e ↑↑ s and e^^s.

Computer Program
The following computer program is written for PowerBasic PBCC and can easily be modified for PBWin, Visual Basic, FORTRAN and any other BASIC-like compiler.

This glossary translates the computer program functions into the mathematical functions of the above analysis when e(x) = exp(x) – 1.

Program	Analysis
DeLog(x)	A^’(x) = a(x) = 1 / Φ(x)
AsympDeLog(x)	1 / asymptotic power series for Φ(x)
eLog(x)	A(x)
eExpInteger(n, x)	e_n(x)
eExp(x)	A _{_– 1}(x)
ExpInteger(n, x)	exp_n(x)
ExpReal(s, x)	exp_s(x)

' The Fourth Operation
' ฉ 2018
'------------------------------------------------------------

'Double precision: plus or minus 4.19E-307 to 1.79*E+308 with 15 or 16 digits of precision
Macro Precision = Double
Macro expRealTol        = .00000001
Macro eExpTolNewton     = .00000001
Macro eExpTolBisection1 = .00000001
Macro eExpTolBisection2 = .00000001
Macro CheckQuit = If WaitKey$ = $Esc Then End

Global e, em1 As Precision
'------------------------------------------------------------

%ht = 8                                       'no more than 8, see Dim coef and Assign coef
Function AsympDeLog(x As Precision) As Precision
 Local i As Long
 Local sum As Precision
 Static coef() As Precision :Dim coef(8)      'static so need only assign once

 If coef(2) = 0 Then                          'first call?
  '                     0  1  2    3      4     5       6        7        8
  Array Assign coef() = 0, 0, 1/2, -1/12, 1/48, -1/180, 11/8640, -1/6720, -11/241920
 End If

 If x < 0 Then
  Print "AsympDeLog: out of bounds" :CheckQuit
  Exit Function
 ElseIf x = 0 Then
  Function = 0
  Exit Function
 End If

 sum = 0
 For i = 2 To %ht                   'terms to x^5 give at least five decimal places if x < .05
  sum = sum + coef(i)*x^i
 Next

 Function = 1/sum
End Function
'------------------------------------------------------------

%maxLogRep = 20
Function DeLog(x As Precision) As Precision
 Dim xx(%maxLogRep) As Precision
 Local d As Precision
 Local i, n As Long

 If     x < 0 Then
  Print "DeLog: out of bounds" :CheckQuit
  Exit Function
 ElseIf x = 0 Then
  Function = 0
  Exit Function
 End If

 i = 0
 xx(0) = x
 Do
  Incr i :If i > %maxLogRep Then Print "DeLog: failed" :Exit Function
  xx(i) = Log(1 + xx(i - 1))
 Loop Until xx(i) < .1                   'anything smaller doesn't affect output before 6 decimal places
 n = i

 d = 1
 For i = 0 To n - 1
  d  = d * (1 + xx(i))
 Next

 Function = AsympDeLog(xx(n)) / d
End Function
'------------------------------------------------------------

Function eLog(ByVal x As Precision) As Precision   'integrate DeLog(t) from 1 to x
 Local cnt As Long                                 'use Romburg's Method of integration
 Local i, k, nmax As Long
 Local sum, h, dx As Precision

 If x < 0 Then
  Print "eLog: out of bounds" :CheckQuit
  Exit Function
 End If

 If x = 1 Then
  Function = 0
  Exit Function
 End If

 cnt = 0                               'shorten interval, at best  x < e - 1,  0 <= eLog(x)
 If x > 0 Then
  Do Until x < em1
   x = Log(1 + x)
   Incr cnt
  Loop
 End If

 dx = .001                             'want nmax so  abs(x - 1)/2^nmax  <=  dx
 nmax = Log2(Abs(x - 1)/dx)
 If nmax = 0 Then nmax = 1
 Dim R(nmax,nmax) As Precision         'typically nmax is no more than 9
                                       'done with dx, all we needed was nmax
 h = (x - 1) / 2
 R(0,0) = h * (DeLog(1) + DeLog(x))

 For i = 1 To nmax
  h = (x - 1) / 2^i
  sum = 0
  For k = 1 To 2^(i-1)
   sum = sum + DeLog(1 + (2*k-1)*h)
  Next
  R(i,0) = R(i-1,0)/2 + h*sum
  For k = 1 To i
   R(i,k) = R(i,k-1) + (R(i,k-1) - R(i-1,k-1))/(4^k - 1)
  Next
 Next

 Function = R(nmax,nmax) + cnt
End Function
'------------------------------------------------------------

Function eExpInteger(ByVal n As Long, ByVal x As Double) As Precision     'e sub n (1)
 Local i As Long, a As Precision
 For i = 1 To n
  x = Exp(x) - 1
  If x = 1/0 Then Exit Function                    'Function = 0, flag for overflow
 Next
 Function = x
End Function

Function eExp(ByVal y As Precision) As Precision   'inverse of eLog, y >= 0
 Local x, xlast, x1, x2 As Precision
 Local n, cnt As Long

 cnt = 0                     'get fractional part of y
 Do Until y < 1
  Decr y
  Incr cnt
 Loop
 Do Until y >= 0
  Incr y
  Decr cnt
 Loop
 'here  0 <= y < 1,  cnt = integral part of original y

 'Try Newton's method first because it usually converges and is fast.
 'It uses the derivative of eLog.
 x = 1
 n = 0
 Do
  xlast = x
  x = xlast + (y - eLog(xlast)) / DeLog(xlast)
  Incr n :If n > 30 Then GoTo BisectionMethod        'hung, use bisection instead
 Loop Until Abs(x - xlast) < eExpTolNewton
 Function = eExpInteger(cnt, (x + xlast) / 2)
 Exit Function

BisectionMethod:
 'bisection method (slow and must specify range of solution)
 x1 = 0              'bracket solution
 x2 = eExpInteger(CLng(y)+1,1)
 Do
  x = (x1 + x2) / 2
  If Abs(eLog(x) - y) < eExpTolBisection1 Or Abs(x1 - x2) < eExpTolBisection2 Then Exit Do
  If Sgn(eLog(x) - y) = Sgn(eLog(x2) - y) Then
   x2 = x
  Else
   x1 = x
  End If
 Loop
 Function = eExpInteger(cnt, x)
End Function
'------------------------------------------------------------

' n must not exceed 3 or 4 (depends on x) to avoid overflow.
' If n = -1, x must > 0,  if n = -2, x must > 1,  if = -3, x must > e, ...
' in general if n < 0, x must > ExpInteger(-2-n,1)
Function ExpInteger(n As Long, ByVal x As Precision) As Precision
 Local i As Long
 If n >= 0 Then
  For i = 1 To n
   x = Exp(x)
  Next
 Else
  For i = 1 To -n
   If x <= 0 Then
    Exit Function            'Function = 0, flag for overflow
   Else
    x = Log(x)
   End If
  Next
 End If
 Function = x
End Function
'------------------------------------------------------------

Function ExpReal(ByVal s As Precision, ByVal x As Precision) As Precision
 Local p, y, z, zlast As Precision
 Local cnt, n As Long

 cnt = 0                     'get fractional part of s
 Do Until s < 1
  Decr s
  Incr cnt
 Loop
 Do Until s >= 0
  Incr s
  Decr cnt
 Loop
 'here  0 <= s < 1,  cnt = integral part of original s

 If s = 0 Then
  Function = ExpInteger(cnt, x)
  Exit Function
 End If

 If x < 0 Then                'get 0 <= x < e
  x = Exp(x)
  Decr cnt
 End If
 Do Until x < e
  x = Log(x)
  Incr cnt
 Loop

 n = 0
 Do
  zlast = z
  p = eLog(ExpInteger(n, x))
  y = eExp(p + s)
  z = ExpInteger(-n, y)
  Incr n
 Loop Until Abs(z - zlast) < expRealTol Or z = 0

 If z <> 0 Then
  Function = ExpInteger(cnt, z)
 Else
  Function = ExpInteger(cnt, zlast)
 End If
End Function
'------------------------------------------------------------

Function PBMain
 Local w, h, a, b As Long
 Local s As Precision
 Desktop Get Size To w, h
 Console Get Size To a, b
 Console Set Loc (w - a)\2, (h - b)\2
 Console Set Screen 120, 100

 e = Exp(1)                          'global 2.71828...
 em1 = e - 1

 Print "    "; " s", " e sub s"
 Print "    "; "----------------------"
 For s = -1 To 3.5 Step .25
  Print "    "; Format$(s,"* .000;-0.000"), Format$(ExpReal(s, 1)," #.00000")
 Next

 WaitKey$
End Function
'------------------------------------------------------------

Bases Other Than e
What if we wanted to know 2_ฝ instead of e_ฝ? And in general b_s for other bases b > 1. We will try to proceed as before when b equaled e and find functional iterates of b^x.

We will use the same labels for the various functions as before. Note that letting

e(x) = b^x – 1

will not work because in general it has two fixed points, one at zero and another either greater or less than zero depending on whether b is less or greater than e. (When b = e the graph of y = exp(x) – 1 is tangent to the graph of y = x, otherwise they intersect in two points one of which is zero.) The trouble with having more than one fixed point is that then there could be no Abel’s function for the function. Recall that the Abel function A(x) must be defined for x > 0 and

A(e(x)) = A(x) + 1

If for example, e(x) = 2^x – 1, which has a fixed point at zero and at one, then A(x) could not be defined at x = 1. The Abel function will always be singular at a fixed point.

Instead of translating the graph of y = b^x vertically by minus one we need to translate it by an amount making it tangent to y = x. In general, unlike the case when b = e, the tangent point won’t be at the origin. To simplify expressions in what follows let

r = log b

where as usual log is the natural logarithm. Note that

b^x = e^rx

When b = e, r = 1.

Setting the derivative of y = b^x to 1 (the slope of y = x) and solving for x, gives x = – (log r) / r. Again to simplify future expressions let

u = – (log r) / r

Note that

b^u = 1/r

So the point (u, 1/r) on the graph of y = b^x is where the slope is 1. Thus if we let

e(x) = b^x – (1/r – u)

then e(u) = u and that is the only fixed point. The derivative is

e^’(x) = r b^x

Let A(x) be an Abel function for the new e(x)

A(e(x)) = A(x) + 1

and let

a(x) = A^’(x)

be the corresponding derived Abel function. Schr๖der’s equation for the new e(x) is

Φ(e(x)) = r b^x Φ(x)

As before we want to determine the Taylor series of Φ(x) but this time at x = u.

Φ(x) = c_₀ + c_₁ (x – u) + c_₂ (x – u)² + c_₃ (x – u)³ + ... as x → u⁺

Taking the first and second derivatives of each side of Schr๖der’s equation and evaluating at u we can, as before, eventually conclude that Φ^’(u) = 0 and Φ^’’(u) = 0, so that both c₀ and c₁ in the Taylor expansion are 0.

Φ(x) = c_₂ (x – u)² + c_₃ (x – u)³ + ... as x → u⁺

Again we need another coefficient and turn to finding the Taylor expansion of e_s(x), first in powers of (x – u) and then in powers of s.

More or less as before we take the first derivative with respect to x of each side of e_s(x) = e(e_s–1(x)) and evaluate at u to eventually obtain

e^’_s (u) = e^’_s–1 (u)

and so

e^’_s (u) = e^’_t (u)

where t is the fractional part of s. Thus as before e^’_s(u) = 1 when s is a whole number and we impose the condition that this holds for all s.

Taking the derivative with respect to x again and evaluating at x = u eventually leads to

e^’’_s (u) = r + e^’’_s–1 (u)

so that

e^’’_s (u) = n r + e^’’_t (u)

where n and t are the integral and fractional parts of s respectively. Thus when s is a whole number e^’’_s(u) = s r and we impose the condition that this holds for all s.

Thus, with our conditions the Taylor expansion begins:

e_s(x) = (x – u) + 1/2 s r (x – u)² + ... as x → u

Determining the Taylor expansion of e_s(x) viewed as a function of s goes exactly as before and we have

e_s(x) = x + Φ(x) s + c_₂ s² + c_₃ s³ + ... as s → 0

We are trying to find Φ(x) . We proved that its Taylor expansion begins

Φ(x) = c_₂ (x – u) ² + c_₃ (x – u) ³ + c_₄ (x – u) ⁴ + ... as x → u⁺

and if, as before, we can determine just c_₂ we can determine the others. Comparing the two series for e_s(x), as when b = e (and r = 1) we impose the following condition, which supersedes our previous conditions because alone it will allow us to determine Φ(x).

Condition
The second order coefficient in the power expansion of Φ(x) is ฝ r.

Thus we have

Φ(x) = 1/2 r (x – u) ² + c_₃ (x – u) ³ + c_₄ (x – u) ⁴ + ... as x → u⁺

Now we use Schr๖der’s equation to determine the other coefficients. Schr๖der’s equation is

Φ(e(x)) = r b^x Φ(x)

Φ(e^rx – 1/r + u) = r e^rx Φ(x)

Substituting in the power series as we did before when b = e we can deduce that

c₃ = (1/48 u² r⁴ – 1/6 u r³ + 7/24 r² – 1/3 r³)
/ (27/8 r² – 2 r – 7/8 +1/24 u³ r³ + 1/2 u² r² + 3/2 u r)

and

c₄ = (25u⁴r⁹ – 16u³r⁹ – 162u²r⁹ – 256u³r⁸ – 576u²r⁸ + 1134ur⁸ + 1296r⁸ + 630u²r⁷ – 480ur⁷ – 1620r⁷ –294ur⁶ + 48r⁶ + 336r⁵)
/ ( –2u⁷r⁹ + 16u⁶r⁸ – 162u⁴r⁸ + 2048u³r⁸ + 168u⁵r⁷ + 96u⁴r⁷ + 24576u²r⁷ + 165888r⁷ –918u⁴r⁶ – 1536u³r⁶ – 42768u²r⁶ + 73728ur⁶ – 255768r⁶ –3968u³r⁵ + 16128u²r⁵ – 31104ur⁵ + 81408r⁵ + 19344u²r⁴ – 9216ur⁴ + 2304r⁴ – 19008ur³ + 3840r³ + 5208r² + 1296u³r⁷)

and succeeding coefficients can be computed likewise. If b = e then r = 1 and u = 0; and c₃ reduces to – 1/12 and c₄ to 1/48, the values we found before. When b = 2, c₃ = – .053891092257005 and c₄ = .00280571012920368.

Repeated application of

e_-1(x) = logarithm to base b of (x– u + 1/r) = (1/r) log (x – u + 1/r)

will reduce any x > u to as close to u as desired. Given x₀ > u, let x_i = e_-i(x₀), then using Schr๖der’s equation as before we can eventually conclude

Φ(x₀) = rⁿ (x₀ + v) (x₁ + v) ... (x_n-1 + v) Φ(x_n)

where v = 1/r – u. Since u might exceed zero we also need to consider the case x₀ < u. Repeated application of e(x) will increase x₀ to as close to u as desired. Let x_i = e_i(x₀), then

Φ(x₀) = Φ(x_n) / (rⁿ(x₁ + v)(x₂ + v) ... (x_n + v))

Thus we have Φ(x) and consequently A(x) and the continuous iterates of e(x). We can get from continuous iterates of e(x) to those of b^x by the same method as before except that exp and log are replaced with exp and log to the base b.

Afterword
A case can be made that there is no series of operations as presented in the introduction, not when the second operand is fractional. In this view there are only two binary operations, addition and multiplication. Exponentiation is not a binary operation but a unary one (the exponent) with a parameter (the base).

Exponentiation cannot be a true binary operation because it cannot be defined when the base is negative and the exponent varies continuously. As a unary operation its purpose is to turn multiplication into addition; it is a go-between for two binary operations but not a binary operation itself.

If you try to think of exponentiation as a binary operation, its non-associatively is a major drawback. (Its noncommutativity is minor in comparison.)

Another reason exponentiation stands out from addition and multiplication is that it has no geometric construction that is continuous in the exponent. Addition can be thought of as laying two line segments end to end in a straight line, multiplication as the area indicated when line segments are joined perpendicularly. Though x, x squared, x cubed, etc. make sense geometrically – a line segment, square, cube, and so forth in higher dimensions – there is no continuum of dimensions to make sense of x raised to a fraction. (The fractional and continuous geometries of Menger, Birkoff, and von Newman aren’t applicable here because their dimension ranges between zero and one, and also they are restricted to projective geometries.)

The problem with negative bases in exponentiation affects what we have been calling the fourth operation, for there, even when the second operand is integral, each number in the “exponential tower” except the highest is an exponential base.

Some simple fractional (not irrational) exponents are useful in physics, where the denominator of the fraction indicates a root. And in mathematics, when trying to describe the growth of a function as a variable goes to infinity a real (even transcendental) exponent might be useful. Szekeres was interested in continuously iterating e(x) because it would interpolate G. H. Hardy’s logarithmico-exponential scale.