Tetration
The following code, regardless of CAS or programming language, will follow some structure, and will always have the following functions:
|
|
Finds an n-th degree polynomial approximate solution to
a given functional equation.
Usage: FSolve[equation, f[x], {x, x0, n}, {coeff}] |
|
|
Prepares a list of power series coefficients from
an n-th degree Abel linearization of a function.
(Formerly known as prepare.) This function is equivalent
to FSolve[a[f[x]] == a[x] + 1, a[x], {x, x0, n}, {y0}]
Usage: Linearize[f[x], {x, x0, n}] |
|
|
Finds the super-logarithm base x of z using the coefficients. |
|
|
Finds the y-th tetration of x using inverse function methods. |
Recently I have been using Mathematica more than Maple, so the most recent version of the code is available for Mathematica only. If I have time or get funding I may take the time to port it over to Maple someday. Until then I would recommend using my most recent package with a trial version of Mathematica if you cannot afford the full version.
(*
** 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]]]]
(* 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];
## 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;
## 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):