<< Computer AlgorithmsUse  Ctrl +  or  –  to enlarge or reduce text size.

The Fourth Operation
- tetration from a real perspective -

Introduction  · Groping for a Solution  ·  Uniqueness  ·  A Related Problem  ·  Facts About Functions  ·  Asymptotic Power Series  ·  To Proceed  ·  Uniqueness Again  ·  The Original Problem  ·  Computer Program  ·  Is it Unique?  ·  Bases Other Than e  ·  Review and Complex Numbers  ·  Other Terminology and Notation  ·  Afterword

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:
AB    =     A · A · A · ...
       B times

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

AB   =		    AAA ···
      B times

We chose the grouping so that the exponentiations are done from top to bottom,  A to the (A to the (A to the ...)).  Any other grouping would collapse, at least in spots, reverting to the third operation, because  (X to the X) to the Y  equals  X to the X·Y.

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:
A1  =  A
AB + 1  =  A · AB
The fourth operation:
A1  =  A
AB + 1  =  AAB
The  nth  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 two operations, addition and multiplication, to include all real numbers. The third, exponentiation, can be substantially so extended. The only restriction is that the base,  A  in  AB,  be nonnegative; this is necessary to ensure that the operation is well-defined and continuous as the exponent,  B,  varies.

The problem we investigate here is how to define the fourth operation,  AB  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.

For now, it will make some formulas simpler if we consider just the case  A = e,  the base of the natural logarithms.
e0  =  1
e1  =  e  =  2.71828…
e2  = ee  =  15.15426...
e3  = eee  =  3814279.10476…
etc.
What is
es
for any real number  s ?

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’ point of view differs from ours in that his primary interest is (see below for definitions) in  exps(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  exps(1)  where  s  is the variable. Jumping ahead, note that Szekeres claims to have proved that the second order coefficient in the power series 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)  =  ex
and place a subscript on  exp  to denote the level of its iteration. The functional iterates of  exp(x)  are closely related to  es.  First observe that
exp0(x)  =  x
exp1(x)  =  exp(x)
exp2(x)  =  exp(exp(x))
exp3(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.  Setting  x = 1  in  expn(x)  gives us the fourth operation without the “socket”  x:
expn(1)  =  ee...e, n times  =  en
 
Thus our goal is to define the function
exps(x)
for all real numbers  s,  such that
exp 1 (x)  =  exp (x)
exp s (exp t (x) )  =  exp s+t (x)
for real  s  and  t.  Then
exps(1)
will answer our original question, what is  es,  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    mul2(x)  =  e2 x
mul (mul (mul (x)))  =  e e e x    or    mul3(x)  =  e3 x
etc.  In this case we know the continuous iterate of mul, it is
muls(x)  =  es x
so that
muls(1)  =  es
Iterates of the function  mul(x)  where  x  is a sort of socket for the first (viewed from the right) 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
AB  =  C
one operation is required to obtain the base  A  and another the exponent  B.  Taking a root gives the base
A  =  Bth root of C  =  B√ C
and taking a logarithm gives the exponent:
B  =  logarithm (base A) of C  =  logAC
where contrary to usage later the subscript  A  here 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
C1/n  =  n√ C
and
Cm/n  =  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  ri  that converges to  s  as  i → ∞  (here the subscript on  r  is an index, not the fourth operation):
lim  ri  =  s
 i → ∞
and defining
Cs  =  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.

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)  =  ∫1x 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 (xp)  =  p f (x)
And then for any whole number  q
q f (x1/q)  =  f (x(1/q)q)  =  f (x)
so that
f (x1/q)  =  1/q f (x)
Thus for any rational number  r
f (xr )  =  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).  We choose the symbol  e  advisedly because it is Euler’s constant.  Then
f (er )  =  r
Thus  f (x)  is the logarithm for exponentiation using the base  e  and we have a very interesting formula for it:

  log x  =  ∫1x 1/ t  dt  

And we have  ex  as well because it is the inverse of  log x.  (We write  log  for the natural logarithm function which in Calculus textbooks is usually written  ln .)

Note that the derivative of  log x  is:
log’x  =  1/x
a simple function independent of  log x  and of  ex.  Thus we can use the function  1/x  to determine the function  log x.  This contrasts with the derivative of  ex , which is  ex  again.

Though the second procedure (using the logarithm function, 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 just 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  es  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
es  =  any function of  s  that is one when  s  is zero
then use the fact that
es+1  =  ees
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
es  =  expn(et)
What is wanted is a condition on our definition of  es  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 properties that recommend it:  (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  ex  that is natural and arguably the best is a difficult problem to attack directly, so instead at first we will try to find such an iterate for a slightly different function
e(x)  =  ex – 1
This function has a fixed point, the only one, at zero:
e(0)  =  0
It is much easier to continuously iterate a function with exactly one 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).  Now please stay out of the Algorithms section.

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

Finding the “logarithm” for  es(1)  is simpler than finding  es(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(es)  =  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 for  f (x),  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(f2(x))  =  A( f ( f (x) ) )  =  A( f (x) ) + 1  =  A(x) + 2
and in general
A( fn(x) )  =  A(x) + n
Thus if we had such a function  A  we could define  fs  for real  s  as the solution to
A( fs(x) )  =  A(x) + s
that is,
fs (x)  =  A -1( A(x) + s )
where  A -1  is the inverse of  A.  Such an  fs (x)  would have the properties
f0 (x)  =  x
fs (fr (x))  =  fs + r (x)
Note that
fs (1)  =  A -1(s)

That works in reverse. If  fs(x)   is a continuum of iterates for  f (x),  that is,  f0(x) = x  and  fs (fr (x)) = fs + r (x) , then  fs  can be obtained from the inverse of the function
B(s)  =  fs(1)
as above.  First replace  s  with  s + r :
B(s + r)  =  fs + r (1)  =  fs(fr(1))
Given  x  let  r = B -1(x).  Then
B(B -1(x) + s)  =  fs(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)  =  ∫1x 1 / Φ(t)  dt  

which generalizes the famous formula for the natural logarithm (third operation).  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) = ex – 1 .

If we take the derivative of each side of Abel’s equation, after some manipulation we get
Φ(e(x))  =  f ’(x) Φ(x)
This is called 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
fs(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)  x2  +  1/6 f ’’’(x) x3  + ...
as  x → 0
We will explain what asymptotic means in the next section. The general term is
1 / i ! f (i)(0)  xi
where  f (i)(x)  denotes the  ith  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 and may or may not converge.

Asymptotic Power Series

Suppose we are given, out of the blue, an infinite series
a0  +  a1(x – x0)  +  a2(x – x0)2  +  a3(x – x0)3  +  ...
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
Sn(x)  =  a0  +  a1(x – x0)  +  a2(x – x0)2  +  ...  an(x – x0)n
then
lim Sn(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 alone could be used to define the function it converges to.

The series always converges at  x = x0,  but except for that value it might not converge. If it doesn’t converge anywhere besides at  x0  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  x0,  any number of its derivatives exist there. Then we can form the Taylor expansion of  F  at  x0
F(x0)  +  F’(x0) (x – x0)  +  1/2 F’’(x0) (x – x0)2  +  ...  1 / i ! F(i)(x0)(x – x0)i  +  ...
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
Sn(x)  =  F(x0)  +  F’(x0) (x – x0)  +  1/2 F’’(x0) (x – x0)2  +  ...  1 / i ! F(n)(x0)(x – x0)n
The sequence might diverge, that is, if we fix  x ≠ x0  and let  n → ∞  the limit of  Sn(x)  might not exist. But even so the series tells something about  F  near  x0. Instead of fixing  x  and letting  n  get arbitrarily large, if we fix  n  and let  x  get arbitrarily close to  x0,  then by definition of the derivative  Sn(x)  will approximate  F(x).

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

If the property allows the latter to approach  x0  arbitrarily closely then the asymptotic series, despite the fact that it doesn’t converge, can be used to compute  F(x)  for any  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 → x0.  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  x0  determines its value far from  x0.

Again, convergence is concerned with  Sn(x)  for fixed  x  as  n → ∞  whereas being asymptotic is concerned with  Sn(x)  for fixed  n  as  x → x0  when the series is the Taylor expansion of a function that has other properties besides Taylor coefficients. Convergence depends on only the coefficients, asymptoticity 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)  =  ex – 1
It has a fixed point at  0.  Its power series expansion about  0  is
e (x)  =  x  +  1/2 x2  +  1/6 x3  + ...
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) = ex
Φ(e(x))  =  ex Φ(x)

We will try to find the Taylor series for  Φ(x)  at  0.  We can expect only that it converges asymptotically. Then, as we shall see later, knowing  Φ(x)  when  x  is near  0  will allow us to compute it for any  x .
Φ(x)  =  c0  +  c1 x  +  c2 x2  +  c3 x3  +  ...      as  x → 0+
It is fairly easy to show that both  c0  and  c1  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)  =  c2 x2  +  c3 x3  +  c4 x4  +  ...      as  x → 0+
where we have yet to determine the remaining coefficients   c2,  c3,  c4,  ...

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  c2,  c3,  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  ci  in order to solve for the others. If we knew  c2  we could solve for c3, then for  c4,  etc.  We will explain in detail when we have figured out what  c2  should be and actually perform the procedure.

Finding  c2  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  es(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 determine  c2  in  Φ(x).  Indeed we can determine the entire series but it will be easier to find it using the above method.

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

First consider  es(x)  as a function of  x  with  s  as a parameter.  Consider its Taylor expansion at  x = 0
es(x)  =  c0  +  c1 x  +  c2 x2  +  c3 x3  +  ...      as  x → 0
The coefficients  ci  are functions of  s ,  namely the  ith derivative of es(x)  with respect to  x  divided by  i !  evaluated at  x = 0 .

Since  e0(x)  =  x  for all  x,  including when  x = 0,  it must be that  c0 = 0  and we can write
es(x)  =  c1 x  +  c2 x2  +  c3 x3  +  ...
Now  c1  is the first derivative of  es(x)  evaluated at  x = 0.  To try to find the first derivative we use a trick:  before taking the derivative break  es(x)  into  e(es–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.
es(x)  =  e(es–1(x))
e’s (x)  =  e’(es–1(x))  e’s–1 (x)
Then since  es–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  es(x)  that generalizes that formula:
Condition #1
e’s (0)  =  1
With that condition  c1 = 1  independent of  s.  Thus the power series for  es(x)  begins:
es(x)  =  x  +  c2 x2  +  c3 x3  +  ...

Now try the same trick with the second derivative. If we can find it  c2  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  D2  as well as double apostrophe denote differentiation twice.
D2[ es(x) ]   =  D[ e’(es–1 (x) ) e’s–1 (x) ]
=  e’’(es–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  es–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  es(x)  that generalizes that formula:
Condition #2
e’’s (0)  =  s

Thus we have the power series for  es(x) beginning:
es(x)  =  x  +  1/2 s x2  +  c3 x3  +  ...      as  x → 0+
We will not need any more conditions. The rest of the coefficients can be determined by breaking  es+1(x)  two ways
e(es(x))  =  es(e(x))
and plugging in the known power series for  e(x)  and what we know of the power series for  es(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. We already have
c1  =  1
c2  =  1/2 s
The next three coefficients are
c3  =  – 1/12 s  +  1/4 s2
c4  =  1/48 s  –  5/48 s2  +  1/8 s3
c5  =  –  1/180 s  +  1/24 s2  –  13/144 s3  +  1/16 s4
so that
es(x)**  =          x  +  1/2 s x2  +  (– 1/12 s  +  1/4 s2) x3  +  (1/48 s  –  5/48 s2  +  1/8 s3) x4  +  (–  1/180 s  +  1/24 s2  –  13/144 s3  +  1/16 s4) x5  +  ...
as  x → 0
We could use this to compute  es(x)  but determining enough terms of the series for accuracy would be difficult. We will continue to focus on finding  Φ(x).  (See the beginning of this section “To Proceed.”)

For any finite number of terms we can rearrange the above series in powers of  s.
es(x)**  =              x  +  (1/2 x2 – 1/12 x3 + 1/48 x4 – 1/180 x5 + ...) s  +  (1/4 x3 – 5/48 x4 + 1/24 x5 + ...) s2  +  (1/8 x4 – 13/144 x5 + ...) s3  +  (1/16 x5 + ...) s4 + ...
still as  x  goes to  0,  not  s .

Now consider  es(x)  as a function of  s  with  x  as a parameter.  Let its Taylor expansion at  s = 0  be
es(x)  =  c0  +  c1 s  +  c2 s2  +  c3 s3  +  ...      as  s → 0
The coefficients  ci  are functions of  x ,  namely the  ith derivative of es(x) with respect to  s  evaluated at  s = 0  divided by  i ! .

Since  e0(x) = x  right away we know  c0 = x .

Now  c1  is the derivative of  es(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
es(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.
es’ (x)  =  1 / A’(A-1(A(x) + s) )
=  1 / A’(es(x))
=  1 / a(es(x) )
Thus since  e0(x) = x
e0’ (x)  =  1 / a(x)  =  Φ(x)
So the power series in  s  begins
es(x)  =  x  +  Φ(x) s  +  c2 s2  +  c3 s3  +  ...      as  s → 0
It’s easy to get all the higher derivatives with respect to  s  from the formula for the first derivative,  es’ (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 the fact) the exact formula for the  ith  derivative of  es(x)  with respect to  s  when  i > 1 is
– (–1)i (2i–3)!! Φ2i–1(es(x))
where !! indicates the double factorial:  N!!  =  N (N – 2) (N – 4) (N – 6) ... where the last factor is  1  or  2  depending on whether  N  is odd or even; in the case here it’s always odd. The  ith  above derivative divided by  i!  and evaluated at  s = 0  gives  ci  for  i > 1
ci  =  – (–1)i (2i–3)!! / i! Φ2i–1(x)

We are trying to find  Φ(x) .  We proved that its Taylor expansion begins
Φ(x)  =  c2 x2  +  c3 x3  +  c4 x4  +  ...      as  x → 0+
and if we can determine just  c2  we can determine the others.  Now we found series for  es(x)  in two ways:  in powers of x and in powers of s, and the second involved  Φ(x).
es(x)  =  x  +  1/2 s x2  +  [higher powers of s and x]            as  x → 0
es(x)  =  x  +  Φ(x) s    +  [higher powers of Φ(x) and s]      as  s → 0
The first can be written (see above)
es(x)  =  x  +  (1/2 x2 – 1/12 x3 + 1/48 x4 + ...) s  +  [higher powers of s and x]            as  x → 0
for  s < 1  and so as  s → 0.  Equating coefficients of the first power of  s  we conclude that
Φ(x)  =  1/2 x2  –  1/12 x3  +  1/48 x4  +  ...      as  x → 0+
We obtained this series by computing a series in two variables. There is a direct way to find the coefficients, as we pointed out at the beginning, if we know just the second order one, and now we know it is  1/2.

Consider Schröder’s equation
Φ(e(x))  =  ex Φ(x)
The power series for  ex  is
ex  =  1  +  x  +  1/2 x2  +  1/6 x3  + ...
and for  e(x)
e (x)  =  x  +  1/2 x2  +  1/6 x3  + ...
Begin by thinking of  Φ(x) as
Φ(x)  =  1/2 x2  +  c x3  + [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 x2  +  (c + 1/2) x3  +  (3/2 c + 7/24) x4  +  [higher powers of x]
and the right
1/2 x2  +  (c + 1/2) x3  +  (c + 1/4) x4  +  [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 x2  –   1/12 x3  +  c x4  + [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  x5
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 x2  – 1/12 x3  + 1/48 x4  – 1/180 x5  + 11/8640 x6  – 1/6720 x7  ...   
  as  x → 0+ 

again obtaining the series for the coefficient of the first degree of  s  in the  x  power series of  es(x).  The fact that the second order coefficient in  Φ(x)  is the same as the corresponding coefficient in  es(x)  entails that all of them are the same.  Expanding  Φ(x)  directly is much easier than expanding  es(x);  it has the advantage of being a function of only one variable.

The coefficients continue (order 8 to 13):  11 / 241920,  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 → ∞
Thus suppose we are given  x = x0 > 0.  Here a subscript on  x  is a sequential index, not the fourth operation. Then the sequence
x0
x1  =  e–1(x0)
x2  =  e–2(x0)
...
xn  =  e– n(x0)
can be continued far enough  (n  can be made large enough)  so that  xn  is as small as desired.

Now apply Schröder’s equation
Φ(e(x))  =  ex Φ(x)
to x1
Φ(x0)  =  Φ(e(x1))  =  ex1 Φ(x1)
then to x2
Φ(x1)  =  Φ(e(x2))  =  ex2 Φ(x2)
Thus
Φ(x0)  =  ex1 ex2 Φ(x2)
And so on through xn
Φ(x0)  =  ex1 ex2 ... exn Φ(xn)
Since  xi+1  =  e–1(xi)  =  log(1 + xi)  we have

  Φ(x0)  =  (1 + x0) (1 + x1) ... (1 + xn–1) Φ(xn)  

We can use this equation to, in a manner of speaking, drop any  x = x0 > 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)  =  ∫1x 1 / Φ(t)  dt
Then the functional iterates of  e(x)  are
es(x)  =  A – 1(A(x) + s)
Note that
es(1)  =  A – 1(s)
And because  es(1)  =  en(es – n(1) )  we have
A – 1(s)  =  en(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 “c2 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) / x2  exists and  ≠ 0
 x → 0
Call this constant  c .

Suppose there were 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 = x0  and consider the sequence
yn  =  e– n(x0)
It goes to  0  (see earlier text). We will use Schröder’s equation for each  Φ  to express yn in terms of  x0. 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  Φ
Φ(yn)  =  e–n Φ(x0)
Φ*(yn)  =  e–n Φ*(x0)
Thus their ratio doesn’t depend on  n.  Call the ratio r.
r  =  Φ*(yn) / Φ(yn)  =  Φ*(x0) / Φ(x0)
so that
Φ*(yn)  =  r Φ(yn)
At this point it looks like  r  depends on  x0.

By the condition we know
lim  Φ(yn) / yn2  =  c
 n → ∞
where  c ≠ 0.  Hence
lim  Φ*(yn) / yn2  =  r / c
 n → ∞
The left side is a constant so  r  is the product of two constants. Thus  r  is constant and independent of  x0,  and
Φ*(x) / Φ(x)  =  r
or
Φ*(x)  =  r Φ(x)
for every  x.

Now recall that for each  Φ  and all  x
∫1x 1 / Φ(t)  dt  =  A (x)
Letting  x = e(1)
∫1e(1) 1 / Φ(t)  dt  =  A (e(1))  =  1
Thus
∫1e(1) 1 / Φ*(t)  dt  =  1 / r ∫1e(1) 1 / Φ(t)  dt
or
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.

First some notation.  We want
exp n (exp m (x) )  =  exp n+m (x)
to hold for all integers so let  exp-1  be the inverse of  exp:
exp-1(x)  =  log(x)
exp-2(x)  =  log(log(x))
etc.

We need define the fractional iterates  exps(x)  only for  0 ≤ s < 1  because  exps(x)  for larger values of  s  can then be found by repeated application of
exps+1(x)  =  exp(exps(x) )
and for smaller values by repeated application of
exps–1(x)  =  exp–1(exps(x) )  =  log(exps(x) )
So in groping for a definition we will suppose  0 ≤ s < 1.

We assume (that is, we require of our definition) that
expr(exps(x) )  =  expr+s(x)
for  r  and  s  fractional just as when they are integral, and since the base, Euler’s constant, is greater than one we assume that  exps(x)  is an increasing function of  s.

Recall that  e(x) = exp(x) – 1  and  es(x)  is its factional iteration.  We make the following further assumption on  exps(x)
Assumption
If  0 ≤ s < 1  then
exps(x)  ≥  es(x).

Consider the identities
exp0(x)  =  e0(x)  +  0
and
exp1(x)  =  e1(x)  +  1
The first just says  x = x  and the second  exp(x) = e(x) + 1.  Since  exps(x)  and  es(x)  are both increasing functions of  s, with the above assumption:

If  0 ≤ s < 1  then
exps(x)  =  es(x)  +  ε
where  0 ≤ ε < 1.

Now let  n  be a whole number and consider
exps(x)  =  exp–n(expn+s(x) )  =  exp–n(exps(expn(x) ) )
Since  exp–n  is  log  iterated  n  times,
exps(x)  =  logn(exps(expn(x) ) )
The idea is to replace the  exps()  on the right with  es()  and let the iterated exp / log wash out the difference as  n  grows large.

By our “ε lemma” above, replacing  x  with expn(x)
exps(x)  =  logn(es(expn(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
logn(es(expn(x) ) )
goes to  0  as  n → ∞.

This is a consequence of the increasingly slow growth of  logn(x)  as either  n  or  x  increases. First note that
log’x  =  1 / x
for  x > 0,  and in general
logn’x  =  1 / x  log x  log2 x  log3 x ... logn–1 x
for  x  large enough that the arguments of all the logarithms exceed zero. Define the sequence of numbers
an  =  expn –1(x)
Then
logn’an  =  1 / expn–1 x expn–2 x expn–3 x ... exp0 x   →   1
as  n → ∞.  And since  logn’x  is a decreasing function of  x  the above limit is true for any sequence whose terms are greater than  an,  such as  es(expn(x))  >  es(expn-1(x))  ≥  expn-1(x).

Letting  y = es(expn(x))  we have
0   <   logn(y + ε) – logn(y)   <   ε logn’(y)
and the right side goes to zero as  n  goes to infinity.

Thus
exps(x)  =  lim  logn(es(expn(x) ) )
                  n → ∞
or using ⚬ to indicate functional composition:

  exps(x)  =  lim  logn ⚬ es ⚬ expn(x)  
                    n → ∞

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

exps(x)  thus defined is a continuum of iterates for  exp(x).  Even for finite  n  we have
logn⚬ e0 ⚬ expn  =  logn ⚬ I ⚬ expn  =  logn ⚬ expn  =  I
and considering
logm⚬ er ⚬ expm ⚬ logn⚬ es ⚬ expn
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
logm⚬ er ⚬ es ⚬ expm  =  logm⚬ er+s ⚬ expm
so that in the limit
expr(exps(x) )  =  expr+s(x)
as desired.

It’s interesting to compare  exps(x)  and  es(x)  as  x  increases. The following table shows their difference  (the  ε  above)  to three decimal places when  s  is between  0  and  1.

  \    x  
     \    
  s   \  		12345678910  …    
 0.00 0.0000.0000.0000.0000.0000.0000.0000.0000.0000.000  …
0.250.1730.1140.0640.0310.0130.0060.0020.0010.000+0.000+  …
0.500.3740.2610.1500.0750.0360.0170.0080.0030.0020.001  …
0.750.6330.4980.3400.2170.1340.0800.0470.0270.0150.008  …
1.001.0001.0001.0001.0001.0001.0001.0001.0001.0001.000  …

If  s < 1  the difference goes to zero, and rapidly, with increasing  x .

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. With seven decimal places the sequence stabilizes after  n = 2.

nformula
 0  1.2710274 
 1  1.6100969 
 2  1.6451508 
 3  1.6451508 

As a check that this great web of arithmetic hangs together we compute  exp½(exp½(1) ).  By  n = 2  it has stabilized to seven decimal places at  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         es
 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 the graph of  y  = ex  with the  y axis  shrunk in order to show more of the graph.


Features to note:
ex  is undefined for  x ≤ –2.  (It would involve the logarithm of a non-positive number.)
As  x → –2+,  ex → –∞.
When  x = –1,  ex = 0.
When  x = 0,  ex = 1.
When  x = 1,  ex = e = 2.71828 ... .
There is a point of inflection between  –1  and  0. A graph is fairly straight near a point of inflexion; this graph is fairly straight throughout the entire interval from  x = –1 to 0.  Throughout that interval  ex  differs from  x + 1  by less than  .01 .
Computer Program

The following computer program is written for PBCC and can be modified for 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.

ProgramAnalysis  
 DeLog(x)  A’(x)  =  a(x)  =  1 / Φ(x)
 AsympDeLog(x)    1 / asymptotic power series for Φ(x) 
 eLog(x) A(x)
 eExpInteger(n, x) en(x)
 eExp(x) A – 1(x)
 ExpInteger(n, x) expn(x)
 ExpReal(s, x) exps(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
'------------------------------------------------------------

Is it Unique?

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

Szekeres proved that, assuming that the derived Abel function is well behaved in a certain way near zero (with that assumption the derived Abel function is well-behaved at infinity too), it is unique. (See the 1962 reference for details.) But this is unconvincing as an argument for uniqueness. It might be that assuming a different sort of good behavior would work as well and lead to a different solution.

Bases Other Than e

 What if we wanted to know  2½  instead of  e½?  And in general  bs  for other bases  b > 1.  We will try to proceed as before when  b  equaled  e  and find functional iterates of  bx.

We will use the same labels for the various functions as before. We need to define  e(x)  having a unique fixed point whose percentage difference with  bx  goes to zero (that is, their ratio goes to one)  as  x  goes to infinity. Letting
e(x)  =  bx – 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,  with other  b > 1  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) = 2x – 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 = bx  vertically by minus one we try translating it by an amount that makes 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
bx  =  e rx
when  b = e,  r = 1.

Setting the derivative of  y = bx  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
bu  =  1/r
So the point  (u, 1/r)  on the graph of  y = bx  is where the slope is  1.  Thus if we translate by – (1/r – u),  letting
e(x)  =  bx – (1/r – u)
then  e(u) = u  and that is the only fixed point.  The derivative is
e’(x)  =  r bx

In the special case where the base  b = e1/e,  which is 1.44466786...,  we have  r = 1/e,  u = e  so that  1/r – u = 0,  that is, no translation is necessary,  y = bx  has exactly one fixed point already.  In this case there is no question that the continuous iterates are unique.

(If we abandon for a moment our condition that  b > 1,  then when  0 < b ≤ 1  there is also exactly one fixed point so that again the continuous iterates are unique, and better, the graph of  y = bx  is not tangent to the graph of  y = x  so they are much easier to compute.)

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 bx Φ(x)

As before we want to determine the Taylor series of  Φ(x)  but this time at  x = u.
Φ(x)  =  c0  +  c1 (x – u)  +  c2 (x – u)2  +  c3 (x – u)3  +  ...      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  c0  and  c1  in the Taylor expansion are  0.
Φ(x)  =  c2 (x – u)2  +  c3 (x – u)3  +  ...      as  x → u+
Again we need another coefficient and turn to finding the Taylor expansion of  es(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  es(x) = e(es–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:
es(x)  =  (x – u) + 1/2 s r (x – u)2  +  ...      as  x → u

Determining the Taylor expansion of  es(x)  viewed as a function of  s  goes exactly as before and we have
es(x)  =  x  +  Φ(x) s  +  c2 s2  +  c3 s3  +  ...      as  s → 0

We are trying to find  Φ(x) .  We proved that its Taylor expansion begins
Φ(x)  =  c2 (x – u) 2  +  c3 (x – u) 3  +  c4 (x – u) 4  +  ...      as  x → u+
and if, as before, we know just  c2  we can determine the others.  Comparing the two series for  es(x),  as when  b = e  (and r = 1)  shows that the second order coefficient in the power expansion of  Φ(x)  is  ½ r.
Thus we have
Φ(x)  =  1/2 r (x – u) 2  +  c3 (x – u) 3  +  c4 (x – u) 4  +  ...      as  x → u+

Now we use Schröder’s equation to determine the other coefficients. Schröder’s equation is
Φ(e(x))  =  r bx Φ(x)
or
Φ(erx – 1/r + u)  =  r erx Φ(x)
Substituting in the power series as we did before when  b = e  we can deduce that
c3  =  (1/48 u2 r4 – 1/6 u r3 + 7/24 r2 – 1/3 r3)
          / (27/8 r2 – 2 r – 7/8 +1/24 u3 r3 + 1/2 u2 r2 + 3/2 u r)
and
c4  =  (25u4r9 – 16u3r9 – 162u2r9 – 256u3r8 – 576u2r8 + 1134ur8 + 1296r8 + 630u2r7 – 480ur7 – 1620r7 –294ur6 + 48r6 + 336r5)
          / ( –2u7r9 + 16u6r8 – 162u4r8 + 2048u3r8 + 168u5r7 + 96u4r7 + 24576u2r7 + 165888r7 –918u4r6 – 1536u3r6 – 42768u2r6 + 73728ur6 – 255768r6 –3968u3r5 + 16128u2r5 – 31104ur5 + 81408r5 + 19344u2r4 – 9216ur4 + 2304r4 – 19008ur3 + 3840r3 + 5208r2 + 1296u3r7)
and succeeding coefficients can be computed likewise.  If  b = e  then  r = 1  and  u = 0;  and  c3 reduces to  – 1/12   and c4  to  1/48,  the values we found before.  When  b = 2,  c3 = – .053891092257005  and  c4 =  .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  x0 > u,  let xi = e-i(x0),  then using Schröder’s equation as before we can eventually conclude
Φ(x0)  =  rn (x0 + v) (x1 + v) ... (xn-1 + v) Φ(xn)
where  v = 1/r – u.  Since  u  might exceed zero we also need to consider the case  x0 < u.  Repeated application of  e(x)  will increase  x0  to as close to  u  as desired.  Let  xi = ei(x0),  then
Φ(x0)  =  Φ(xn) / rn(x1 + v)(x2 + v) ... (xn + 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  bx  by the same method as before except that  exp  and  log  are replaced with  exp  and  log  to the base  b.  ( In the special case  b = e1/e  this last step is unnecessary because then  e(x) = b(x). )

Review and a Glance at the Complex Plane

We sought to define  e  raised to itself s times where  s  can be fractional. We added a sort of “socket” variable to turn the problem into defining continuous iterates of the function  exp(x).

Since finding a unique continuous iterate of a function that doesn’t have a fixed point, like exp(x), is difficult, instead we considered the function  exp(x) – 1,  which has exactly one fixed point. Following George Szekeres (1911 – 2005) and patching up gaps in his work, we found by an intricate argument and some reasonable assumptions a unique solution for the continuous iterates of that function.

Then we used that to find a solution for the continuous iterates of  exp(x).  However its uniqueness is in doubt.

So far all we have used is real analysis. Another way to get around the fact that  exp(x)  has no fixed point is to go into the complex plane to find one, that is, consider  exp(z)  where  z  varies over the complex numbers. The trouble is that exp(z) has more than one fixed point, indeed an infinite number, and that leads to singularities and non-uniqueness. There are ways of dealing with the singularities but to get uniqueness seems to require arbitrary conditions. In any case Hellmuth Kneser (1898 – 1973) used one fixed point to construct an analytic solution, and later other people improved on his work. As for numerical results, Kneser’s solution differs from Szekeres’. For example, while  e½ is 1.64515… per Szekeres, it is 1.64635… per Kneser.

There are many Kneser’s iterations one for each fixed point. Ultimately the choice it arbitrary – one fixed point is like another – so Kneser’s primary iteration is not unique. Szekeres’ iteration has the advantage of being easier to compute and it can be argued that it is the best iteration at infinity.

At any rate, even if we failed with  exp(x),  the continuous iteration of exp(x) – 1,  solved here without arbitrary assumptions, is interesting in its own right.

Szekeres was interested in continuously iterating  e(x)  =  exp(x) – 1  in order to interpolate G. H. Hardy’s logarithmico-exponential scale.  a(x), its derived Abel function, grows more slowly than any finite iteration of  log(x)  and its inverse grows more rapidly than any finite iteration of  exp(x),  thus these two functions are respectively below and above Hardy’s scale.

Other Terminology and Notation

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

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

Afterword

A case could be made that there is no sequence of operations as presented in our introduction, that there are only two binary operations, addition and multiplication; that exponentiation is not a binary operation but a unary one (the exponent) with a parameter (the base). As a unary operation its purpose is to turn addition into multiplication; it is a go-between for two binary operations but not a binary operation itself.

Exponentiation cannot be a true binary operation because it cannot be defined when the base is negative and the exponent varies continuously.

Because of the non-associatively of exponentiation it cannot be a group operation, even if zero or all non-positive real numbers are excluded.

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. We run into a problem with exponentiation.  X,  X squared,  X cubed,  etc.  make sense geometrically – a line segment, square, cube, and so forth in higher dimensions – but there is no continuum of dimensions to make sense of  X  raised to a fraction. What, for example, is between a square and a cube? (The fractional and continuous geometries of Menger, Birkoff, and von Newman are no help because their dimensions range between zero and one, and they are projective geometries rather than Euclidean.)

A further difference – among all three of addition, multiplication, and exponentiation – becomes apparent when using the operations to solve practical problems. In practice many numbers are not pure numbers but quantities, that is, they measure something in a certain unit such as meters or kilograms. We will say a pure number has the “null unit” so that all numbers have units. When adding two numbers both must be in the same unit. “You can’t add apples and oranges” as the saying goes. When multiplying there is no such restriction, each number can have any unit. When neither number is pure the result is in a new unit. For example, two meters times three meters is six square meters, two watts times three hours is six watt-hours.

Multiplication can be considered iterated addition only when one of the numbers – the number measuring the degree of iteration – is pure. Adding five apples to itself three times makes 15 apples. But three apples or three oranges as a count for adding makes no sense.

With exponentiation, the exponent must be pure, and if the exponent can vary continuously the base must be pure as well. In general there is no exponentiation of quantities.  (In physics an exponent might be an expression the individual terms of which have a unit but when evaluated the units cancel.)

Thus it would seem that we have gone on a fool’s errand, though an interesting one. Not only is there no fourth operation, there is not even a third operation. Exponentiation should not be viewed as a third binary operation or as any kind of binary operation. It is a unary operation on pure numbers, or rather a continuum of unary operations parameterized by a base. It is a bridge between the two elementary operations on numbers stripped of their units that converts addition to multiplication and whose inverse converts multiplication to addition.

Some simple fractional 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.