[Math] Algorithm for tetration to work with floating point numbers

algorithmsexponentiationhyperoperationtetration

So far, I've figured out an algorithm for tetration that works.

However, although the variable a can be floating or integer, unfortunately, the variable b must be an integer number.

How can I modify the algorithm so that both a and b can be floating point numbers and the correct answer will be produced?

// Hyperoperation type 4:
public double tetrate(double a, double b)
{
    double total = a;
    for (int i = 1; i < b; i++) total = pow(a, total);
    return total;
}

In an attempt to solve this, I've create my own custom power() function (trying to avoid roots, and log functions), and then successfully generalized it to multiplication. Unfortunately, when I then try to generalize to tetration, numbers go pear shaped.

I would like an algorithm to be precise up to x amount of decimal places, and not an approximation as Wikipedia talks about.

Best Answer

When I studied the various known matrices of combinatorical numbers I also looked at the following idea: what if we discuss functions with a set of results, not only one number? So for instance the concept of sine and cosine gets some special charme if we look not only at f(x) = sin(x), g(x)=cos(x) but at 2x2-matrices containing cos() and sin() and the input of matrix-multiplications has two parameters and the output as well. This helps much to understand algebra with complex numbers, for instance.
After that I thought: what if we look at matrices which allow not only to look at a function $f(x)=a_0 + a_1x + a_2x^2 + a_3x^3 + ...$ involving all the powers of $x$ as input, but also spit out all powers of $f(x)$ in the same manner. Then for instance

$$ V(x) \cdot P = V(x+1) $$ where $V(x)$ means the vector $V(x)=[1,x,x^2,x^3,x^4,...]$ and $P$ means the upper triangular Pascal matrix. Of course $$(V(x) \cdot P )\cdot P = V(x+2)$$ and of course it is easily to see that $$V(x) \cdot P^h = V(x+h)$$ first for integral powers of $P$ which any programmer can implement - and if you download Pari/GP you can just use the matrix-language of it.
It needs a certain ingenious mind to find out how to define fractional powers of $P$ to make the fractional operation of addition possible...

After that one can find, that a very well known combinatorical matrix $B$ performs $$ V(x) \cdot B = V(exp(x))$$
and logically $$ V(x) \cdot B^h = V(exp^{oh}(x))$$
for integer $h$ - but which means implementing tetration to integer iterates $h$.
To find now the fractional-h power of $B$ -which is your question- is not so simple, and using matrices practically means to truncate them to finite size.

Well, anyway meaningful approximations can be done and truncated versions of $B$ can be diagonalized and general powers of it can then be found by simply computing the general powers on each of the scalar numbers in the diagonal - and so one can introduce fractional powers of $B$ in the formula $$ V_{32}(x) \cdot B_{32x32}^h \approx V_{32}(exp^{oh}(x))$$
(for for instance $32$ and $32x32$ sized vectors and matrix) and this implements then the fractional iteration of the exponential, aka tetration.

The matrices $P$ and $B$ of which I talk here are known as "Carleman-matrices" and this whole concept can be applied to fractional iteration of many functions for which you have a power series. (Well, we'll have convergence issues and much more complication and non-uniqueness and what not - but that's not the focus of this answer)
So using "mateigen" in Pari/GP for the diagonalization allows the following Pari/GP code:

n=32
V(x,dim=n) = vector(dim,r,x^{r-1})      \\ define the function for generating V(x)
b = sqrt(2)
bl = log(b)
B = matrix(n,n,r,c,bl^(r-1)*(c-1)^(r-1)/(r-1)!)  \\ define the Carlemanmatrix 
    \\ for   V(x) * B = V(b^x) 
Y = V(1) * B  \\ in Y we'll have approximately V(b^x) 
    \\where the first two entries have the best approximation to the expected
    \\ value
    \\ now diagonalize
M = mateigen(B) \\ needs high precision of say 800 or 1600 digits for all operations
W = M^-1 
D = diag(W*B*M) 
   \\ after this we can do B^h = M * D^h * W 
   \\  and  D^h can be done by d_k,k^h on the elements of the diagonal
dpow(A,h=1)= for (r=1,#A, A[r,r]=A[r,r]^h);return(A) \\ define a function for powers of D      

BPow (h)= M*dpow(D,h)*W  \\ define a function for the fractional power of B

\\ after this we can do
Y = V(1) * BPow(0.5) 
\\ and have in Y[2] a nice approximation for the half-iterate or b^^0.5
\\where b=sqrt(2)

I you want to experiment with this I recommend to use dimension $n=16$ first (with much poorer approximation) but already you'll need default(realprecision,200) or so. To do mateigen on matrix $B$ of size $32x32$ one must compute $b$, $bl$, $B$ and so on already to even $1600$ or $2000$ internal digits precision and the mateigen can need several minutes to complete. I tried this also with size $n=64 $ and this was really a full afternoon computation... but is a nice exercise anyway.

Remark: I've discussed this a bit in the tetrationforum and have made an article which compares a couple of such naive approaches to fractional exponential-towers. The above method I've called there "polynomial approach" because using finite matrices without taking care for compatibility with the theoretically infinite sized matrices and their fractional powers means just involve polynomials for that approximative solutions. See the index and the link to the essay

Related Question