<< Computer AlgorithmsUse  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:
AB    =     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):

AB   =
    AAA ···
     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:
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 three operations to include all real numbers (except that in the case of exponentiation the base,  A  in  AB,  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,  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  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.
e0  =  1
e1  =  e  =  2.71828
e2  = ee  =  15.15426...
e3  = eee  =  3814279.10476

Again, what is
e˝
e  to itself one half times? And in general what is
es
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  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 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)  =  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.  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”:
expn(1)  =  ee...e, n times  =  en
 
Thus our goal is to define the function
exps(x)
for all real numbers  s,  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 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
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  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 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)  =  ∫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)
Also
q f (x1/q) = f (x1/qq) = f (xq1/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). 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.

Note that the derivative of  log x  is:
logx  =  1/x
a simple function independent of  log x  and of  ex.  Thus we can use  1/x  to determine  log x.  This contrasts with the derivative of  ex , which is  ex  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  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 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  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 a continuous iterate for a function not much different
e(x)  =  ex – 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).  es(x)  is a function and the function it iterates is  exp(x) – 1,  yet if we were to write the function without its argument,  es,  it would look like the fourth operation — e^(e^(...^e)...)  s  times — rather than the  sth  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  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, 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
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)
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
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 but might 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 would stand alone and 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 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  x  in terms of its value at points near  x0.

If the property allows us to approach  x0  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 → 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  determine 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 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)  =  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, nonetheless, 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 make a plausible choice for  c2  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  c2 .

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 finding 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))
es (x)  =  e(es–1(x)) es–1 (x)
Then since  es–1(0) = 0  and  e (0) = 1
es (0)  =  es–1 (0)
Repeating that formula within itself on the right we get down to
es (0)  =  et (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  es (0)  in general. However in the particular cases when  s  is a whole number, that is when  t = 0,  we have determined it:
en (0)  =  1
for every whole number  n.  We shall impose the natural, simplest, condition on  es(x)  that generalizes that formula:
Condition #1
es (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) ) es–1 (x) ]
=  e’’(es–1 (x) ) es–1 (x) es–1 (x)  +  es–1 (x) e’’s–1 (x)
Setting  x = 0  and using the fact that  es – 1(0) = 0  and  es – 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. The next three coefficients are
c3  =  1/4 s2 – 1/12 s
c4  =  1/8 s3 – 5/48 s2 + 1/48 s
c5  =  1/16 s4 – 13/144 s3 + 1/24 s2 – 1/180 s
so that
es(x)22  =               x  +  1/2 s x2  +  (1/4 s2 – 1/12 s) x3  +  (1/8 s3 – 5/48 s2 + 1/48 s) x4  +  (1/16 s4 – 13/144 s3 + 1/24 s2 – 1/180 s) x5  +  ...
as  x → 0
For any finite number of terms we can rearrange in powers of  s.
es(x)22  =               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 it) 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) ... 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  ith  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 formulas for  es(x)  in two ways 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  +  ...      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 x2 – 1/12 x3 + 1/48 x4 + ...      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  c2  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 x2  +  c3 x3  +  c4 x4  +  ...      as  x → 0+

Now to find the remaining coefficients.  Consider Schröder’s equation
Φ(e(x))  =  ex Φ(x)
We know the power series for  ex
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)2  =    1/2 x2  1/12 x3  + 1/48 x4  1/180 x5  + 11/8640 x6  1/6720 x7  11/241920 x8  ...
as  x → 0+
Remarkably, this appears to be the first degree coefficient of the power series of  es(x)  in  s.  Positing that the second order term in  Φ(x)  be the same as the corresponding term in  es(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  “c2.”  In any case, 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 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 = 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 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 = 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.

We need define  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 at first suppose  0 ≤ s < 1.  At this point all we impose on  exps(x)  is that it be an increasing function of  s  and of  x,  and
expr(exps(x) )  =  expr+s(x)
for r and s fractional just as when they are integral.

Recall that  e(x) = exp(x) – 1.  Consider the identities
exp0(x)  =  e0(exp0(x) ) + 0
and
exp1(x)  =  e1(exp0(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  the following condition on  exps(x)  is reasonable.
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 assumption 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
logx  =  1 / x
for  x > 0,  and in general
lognx  =  1 / x log x log2x log3x ... logn–1x
for  x  large enough that the arguments of all the logarithms exceed zero. Define the sequence of numbers
an  =  expn –1(x)
Then
lognan  =  1 / expn–1x expn–2x expn–3x ... exp0  →   1
as  n → ∞.  And since  lognx  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 → ∞
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  exps(x)  differs from  es(x)  when  s  is between  0  and  1.  To three decimal places:

   x  = 12345678910
sεεεεεεεεεε
 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

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 .

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 in
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

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.

nformula
 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         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 a graph of y  = ex  with the scale of 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 = 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  ex  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,  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 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  es  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.

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
'------------------------------------------------------------
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. Note that 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,  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) = 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 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
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 let
e(x)  =  bx – (1/r – u)
then  e(u) = u  and that is the only fixed point.  The derivative is
e(x)  =  r bx

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
es (u)  =  es–1 (u)
and so
es (u)  =  et (u)
where  t  is the fractional part of  s.  Thus as before  es(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 can determine just  c2  we can determine the others.  Comparing the two series for  es(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) 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.

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.