Tetration
  Main Identities Code Graphs Paper  
  History Definitions Error Tables Links  

About the Code

Functions


         The following code, regardless of CAS or programming language, will follow some structure, and will always have the following three functions:

New Code


         This page includes code that is very different than the code found in the paper. The code in the paper requires that preparation be done before-hand, whereas the new code allows you to call the tetrate and superlog functions without any preparation, because it does the preparation on the spot. This is very innacurate, but on the good side, it's fast.
         Another new feature is multiple preparation. You can now call prepare (in Mathematica only) with a list describing the numbers of n you wish to prepare. Instead of solving the systems of equations separately, it combines several steps so that it saves time. The way that it does this, is through partial row-reduction. The partial row-reduction method could be the subject of an entire paper, but simply put, it solves a m by m matrix only upto row n (assuming n < m), thus being able to grab solutions from that matrix as if it were an n by n matrix. Then continue by row-reducing the m by m matrix that is already partially row-reduced. This allows you to prepare(n), and prepare(m) at the same time!
         The reason the partial row-reduction method is useful is because the only way to determine the "accuracy" of the values you get with prepare(n) is to compare them to values at a higher preparation, say maybe prepare(2*n). For small values of n, this is trivial, but when comparing values prepared from n=200 to n=400, every microsecond counts.

Mathematica Code

Code from Paper (v1.0)


(*
** Usage:
**      prep = SuperLogPrepare[n, x];
**      SuperLog[prep, z]    -- gives n-th approx. of slog_x(z) 
**      Tetrate[prep, y]     -- gives n-th approx. of x^^y
**
** Copyright 2005 Andrew Robbins
*)

SuperLogPrepare[n_Integer, x_] := {x, LinearSolve[Table[
    k^j/k! - If[j == k, Log[x]^-k, 0], {j, 0, n - 1}, 
    {k, 1, n}], Table[If[k == 1, 1, 0], {k, 1, n}]]}

SuperLog[v_, z_?NumericQ] := Block[{(*SlogCrit*)}, 
    SlogCrit[zc_] := -1 + Sum[v2, k*zc^k/k!, {k, 1, Length[v2]}]; 
    Which[z ² 0, SlogCrit[v1^z] - 1, 0 < z ² 1, SlogCrit[z], z > 1, 
    Block[{i=-1}, SlogCrit[NestWhile[Log[v1, #]&, z, (i++;#>1)&]]+i]]]

Tetrate[v_, y_?NumericQ] := Block[{(*SlogCrit, TetCrit*)},
    SlogCrit[zc_] := -1 + Sum[v2, k*zc^k/k!, {k, 1, Length[v2]}];
    TetCrit[yc_] := FindRoot[SlogCrit[z] == yc, {z, 1}]1, 2; If[y > -1, 
    Nest[Power[v1, #]&, TetCrit[y - Ceiling[y]], Ceiling[y]],
    Nest[Log[v1, #]&, TetCrit[y - Ceiling[y]], -Ceiling[y]]]]

New Code (v1.1)


(* Usage:
** 
** prep = SuperLogPrepare[x, n];    -- Prepare for slog_x(z)_n
** SuperLog[x, z]                   -- Prepare fast (n=10)
** SuperLog[prep, z]                -- Use local prep
** SuperLog[z]                      -- Use global prep
*)
SuperLogPrepare[x_, n_Integer] := {x, LinearSolve[Table[
    k^j/k! - If[j == k, Log[x]^-k, 0], {j, 0, n - 1}, 
    {k, 1, n}], Table[If[k == 1, 1, 0], {k, 1, n}]]}
SuperLogPrepare[x_, {nMax_Integer}] := SuperLogPrepare[x, {1, nMax}]
SuperLogPrepare[x_, {nMin_Integer, nMax_Integer}] := 
    SuperLogPrepare[x, {nMin, nMax, 1}]
SuperLogPrepare[x_, {nMin_Integer, nMax_Integer, nStep_Integer}] := 
    Block[{m}, m = Table[If[k <= nMax, 
        k^j/k! - If[j == k, Log[x]^-k, 0], If[j == 0, 1, 0]],
        {j, 0, nMax - 1}, {k, 1, nMax + 1}];
    Table[m[[Range[i]]] = RowReduce[m[[Range[i]]]];
        {x, m[[Range[i], nMax + 1]]}, 
        {i, nMin, nMax, nStep}]]
SuperLogLinear[x_, z_] := SuperLog[{x, {}}, z, (#2-1)&]
SuperLogCritical[{x_, terms_}, z_] := -1 + 
    Sum[terms[[k]]*z^k/k!, {k, 1, Length[terms]}]
SuperLog[z_?NumericQ] := SuperLog[$SuperLogPrep, z]
SuperLog[{x_?NumericQ, terms_List}, z_?NumericQ] := 
    SuperLog[{x, terms}, z, SuperLogCritical]
SuperLog[{x_?NumericQ, terms_List}, z_?NumericQ, f_] := Which[
    z < 0,      f[{x, terms}, x^z] - 1,
    z == 0,     -1,
    0 < z < 1,  f[{x, terms}, z],
    z == 1,     0,
    z > 1,      Block[{i = -1}, f[{x, terms}, 
        NestWhile[Log[x,#]&,z,(i++;#>1)&]]+i]]
InternalSLO2[x_, z_] := Block[{dec, mtx, y1, y2},
    mtx = Table[Table[If[k == 12, If[j == 0, 1, 0],
        k^j/k! - If[j == k, Log[x]^-k, 0]], {k, 1, 12}], {j, 0, 10}];
    mtx[[Range[10]]] = RowReduce[mtx[[Range[10]]]];
    y1 = SuperLog[{x, Most[mtx][[All, 12]]}, z];
    y2 = SuperLog[{x, RowReduce[mtx][[All, 12]]}, z];
    If[y1 != y2, dec = Floor[-Log[10, Abs[y2 - y1]]];
        Floor[y2*10^dec]/10^dec, y2]];
InternalSLO3 = Compile[{{x, _Real}, {z, _Real}}, InternalSLO2[x, z]];
SuperLog[x_?NumericQ, z_?NumericQ] /; (1.4 <= x <= 16) := InternalSLO3[x, z]

(* Usage:
**
** prep = SuperLogPrepare[x, n];    -- Prepare for tet_x(y)_n
** Tetrate[x, y]                    -- Prepare fast (n=10)
** Tetrate[prep, z]                 -- Use local prep
** Tetrate[z]                       -- Use global prep
*)
TetrateLinear[x_, y_] := Tetrate[{x, {}}, y, (#2+1)&]
TetrateCritical[prep_, y_] :=
    FindRoot[SuperLogCritical[prep, z] == y, {z, 1}][[1, 2]]
Tetrate[y_?NumericQ] := Tetrate[$SuperLogPrep, y]
Tetrate[{x_?NumericQ, terms_List}, y_?NumericQ] := 
    Tetrate[{x, terms}, y, TetrateCritical]
Tetrate[{x_?NumericQ, terms_List}, y_?NumericQ, f_] := If[y > -1, 
    Nest[Power[x, #]&, f[{x, terms}, y - Ceiling[y]], Ceiling[y]],
    Nest[Log[x, #]&, f[{x, terms}, y - Ceiling[y]], -Ceiling[y]]]
InternalTEO2[x_, y_] := Block[{dec, mtx},
    mtx = Table[Table[If[k == 12, If[j == 0, 1, 0],
        k^j/k! - If[j == k, Log[x]^-k, 0]], {k, 1, 12}], {j, 0, 10}];
    mtx[[Range[10]]] = RowReduce[mtx[[Range[10]]]];
    z1 = FindRoot[SuperLog[{x, 
        Most[mtx][[All, 12]]}, z] == y, {z, 1}][[1, 2]];
    z2 = FindRoot[SuperLog[{x, 
        RowReduce[mtx][[All, 12]]}, z] == y, {z, 1}][[1, 2]];
    If[z1 != z2, dec = Floor[-Log[10, Abs[z2 - z1]]];
        Floor[z2*10^dec]/10^dec, z2]]
InternalTEO3 = Compile[{{x, _Real}, {y, _Real}}, InternalTEO2[x, y]];
Tetrate[x_?NumericQ, y_?NumericQ] /; (1.4 <= x) := If[y > -1, 
    Nest[Power[x, #]&, InternalTEO3[x, y - Ceiling[y]], Ceiling[y]],
    Nest[Log[x, #]&, InternalTEO3[x, y - Ceiling[y]], -Ceiling[y]]]

(* Fast Approximations:
**
** SuperRoot[y, z]  == srt_y(z)     (n=10)
** IterExp[x, y, a] == exp_x^y(a)   (n=10)
** IterLog[x, y, a] == log_x^y(a)   (n=10)
*)
SuperRoot[y_, z_] := FindRoot[Tetrate[x, y] == z, {x, 2}][[1, 2]]
IterExp[x_, y_, a_] := Tetrate[x, y + SuperLog[x, a]]
IterLog[x_, y_, a_] := IterExp[x, -y, a]

(* Globals *)
$SuperLogPrep = SuperLogPrepare[E, 10];

Maple Code

Code from Paper (v1.0)


## Usage:
##      prep := superlog_prepare(n, x):
##      superlog(prep, z);  -- gives n-th approx. of slog_x(z)
##      tetrate(prep, y);   -- gives n-th approx. of x^^y
##
## Copyright 2005 Andrew robbins
##

with(linalg):

superlog_prepare := proc(n::integer, x)	[x, linsolve(matrix([seq([seq(		k^j/k!-`if`(j=k,1,0)*log(x)^(-j), 		k=1..n)], j=0..(n - 1))]),		[seq(`if`(k=1,1,0),k=1..n)])];end proc;

superlog := proc(v, z) local slog_crit;    if not (z::numeric) then return 'procname'(args); end if;    slog_crit := proc(zc) -1 + sum(v[2][k]*zc^k/k!,         k=1..(vectdim(v[2]))); end proc;    piecewise(z = -infinity, -2,        z < 0, slog_crit(v[1]^z) - 1, z = 0, -1,        0 < z and z < 1, slog_crit(z), z = 1, 0,        z > 1, (proc() local a, i; a := evalf(z);            for i from 0 while (a > 1) do a := evalf(log[v[1]](a)); end do;            slog_crit(a) + i; end proc)());end proc;

tetrate := proc(v, y) local tet_crit;    if not (y::numeric) then return 'procname'(args); end if;    tet_crit := proc(yc) local slog_crit;        slog_crit := proc(zc) -1 + sum(v[2][k]*zc^k/k!, 		k=1..(vectdim(v[2]))); end proc;        select((proc(a) evalb(Im(a) = 0 and 0 <= Re(a) and Re(a) <= 1) 		end proc), [solve(evalf(slog_crit(z)) = yc, z)])[1];    end proc;    piecewise(y = -2, -infinity,        -2 < y and y < -1, log[v[1]](tet_crit(y+1)),        y = -1, 0, -1 < y and y < 0, tet_crit(y),        y = 0, 1, y > 0, (proc () local a, i;            a := tet_crit(y - ceil(y));            for i from 1 to ceil(y) do a := v[1]^a; end do;            evalf(a); end proc)());end proc;

New Code (v1.1)


## Usage:#### prep := superlog_prepare(x, n):  -- Prepares for slog_x(z)_n## superlog(x, z);                  -- Prepares fast (n=10)## superlog(prep, z);               -- Uses 'prep'## superlog(z);                     -- Uses 'Global_superlog_prep'##with(linalg):superlog_prepare := proc(x, n::integer)    [x, linsolve(matrix([seq([seq(        k^j/k!-`if`(j=k,1,0)*log(x)^(-j),         k=1..n)], j=0..(n - 1))]),        [seq(`if`(k=1,1,0),k=1..n)])];end proc;superlog_linear := proc(x, z)    superlog([x, []], z, (proc(v, z) z - 1; end proc));end proc;superlog_critical := proc(v, z)    -1 + sum(v[2][k]*z^k/k!, k=1..vectdim(v[2])); end proc;superlog := proc() local prep, x, z, f;    if (nops([args]) = 1) then        if not ([args][1]::numeric) then return 'procname'(args); end if;        superlog(Global_superlog_prep, args, superlog_critical);    elif (nops([args]) = 2) then        z := evalf([args][2]);        if not (z::numeric) then return 'procname'(args); end if;        if not ([args][1]::list) then            x := evalf([args][1]);            prep := [x, linsolve(matrix([seq([seq(                evalf(k^j/k!-`if`(j=k,1,0)*log(x)^(-j)),                 k=1..10)], j=0..9)]),                [seq(`if`(k=1,1,0),k=1..10)])];            superlog(prep, z, superlog_critical);        elif ([args][1]::list) then            superlog(args, superlog_critical);        else return 'procname'(args);        end if;    elif (nops([args]) = 3) then        z := [args][2];        f := [args][3];        if not (evalf(z)::numeric) then return 'procname'(args); end if;        if (([args][1])::list) then            prep := [args][1];            x := [args][1][1];            if not (evalf(x)::numeric) then return 'procname'(args); end if;        else return 'procname'(args); end if;        piecewise(            z < 0,              f(prep, x^z) - 1,             z = 0,              -1,            0 < z and z < 1,    f(prep, z),             z = 1,              0,            z > 1,              (proc()                    local a, i; a := evalf(z);                    for i from 0 while (a > 1) do                         a := evalf(log[x](a));                    end do; f(prep, a) + i;                 end proc)())    else return 'procname'(args);    end if;end proc;## Usage:#### prep := superlog_prepare(x, n):  -- Prepares for slog_x(z)_n## tetrate(x, y);                   -- Prepares fast (n=10)## tetrate(prep, y);                -- Uses 'prep'## tetrate(y);                      -- Uses 'Global_superlog_prep'##tetrate_critical := proc(v, y) local slog_critical;    select((proc(a) evalb(Im(a) = 0 and 0 <= Re(a) and Re(a) <= 1)         end proc), [solve(evalf(superlog_critical(v, z)) = y, z)])[1];end proc;tetrate := proc() local prep, x, y, f;    if (nops([args]) = 1) then        if not ([args][1]::numeric) then return 'procname'(args); end if;        tetrate(Global_superlog_prep, args, tetrate_critical);    elif (nops([args]) = 2) then        y := evalf([args][2]);        if not (y::numeric) then return 'procname'(args); end if;        if not ([args][1]::list) then            x := evalf([args][1]);            prep := [x, linsolve(matrix([seq([seq(                evalf(k^j/k!-`if`(j=k,1,0)*log(x)^(-j)),                 k=1..10)], j=0..9)]),                [seq(`if`(k=1,1,0),k=1..10)])];            tetrate(prep, y, tetrate_critical);        elif ([args][1]::list) then            tetrate(args, tetrate_critical);        else return 'procname'(args);        end if;    elif (nops([args]) = 3) then        y := [args][2];        f := [args][3];        if not (evalf(y)::numeric) then return 'procname'(args); end if;        if (([args][1])::list) then            prep := [args][1];            x := [args][1][1];            if not (evalf(x)::numeric) then return 'procname'(args); end if;        else return 'procname'(args); end if;        piecewise(            y = -2,             -infinity,            -2 < y and y < -1,  log[x](tetrate_critical(prep, y + 1)),            y = -1,             0,             -1 < y and y < 0,   tetrate_critical(prep, y),            y = 0,              1,             y = 1,              x,            y > 0,              (proc ()                     local a, i; a := tetrate_critical(prep, y - ceil(y));                    for i from 1 to ceil(y) do a := x^a; end do;                    evalf(a);                end proc)());    else return 'procname'(args);    end if;end proc;Global_superlog_prep := superlog_prepare(exp(1), 10):

Visited: times, last updated: 2006-02-15, by: Andrew Robbins, contact: and_j_rob(at)yahoo(dot)com