+- Tetration Forum (https://tetrationforum.org)
+-- Forum: Tetration and Related Topics (https://tetrationforum.org/forumdisplay.php?fid=1)
+--- Forum: Computation (https://tetrationforum.org/forumdisplay.php?fid=8)
+--- Thread: PARI/GP code to approximate Kneser tetration(base e) without using schroder/abel (/showthread.php?tid=1825)
PARI/GP code to approximate Kneser tetration(base e) without using schroder/abel - RaeesHarris - 07/27/2026
Hello, for the past few days I have succesfully developed an algorithm in PARI/GP to approximate Kneser tetration for Base \(e\) on \([0,1]\) with precision around ~\( 10^{-15} \) (numerically compared to fatou.gp) which is near the Double Machine Epsilon precision. The algorithm doesn't use standard computational objects like Schroder or Abel because I was frustrated with the difficult of numerically implementing those functions in PARI/GP so instead I have developed another method that does not build those computational objects. This is going to be a bit of a long post, but i will try to explain all of the theory behind this
This is my algorithm to approximate Kneser Tetration for Base \(e\) on \([0,1]\):
Code:
\\ Numerical algorithm to calculate an approximation of Kneser tetration on the interval [0,1]
\\ this is V2 which includes taking a pade approximant to make the error near ~10^-15 on the specificied interval
\\ fixed point for f(x) = exp(x) + b
default(seriesprecision, 30);
L(x) = x-lambertw(-exp(x),-1);
\\ faa-di-bruno reccurence series for the superfunction of f(x) = exp(x) + b.
\\ Works for b < -1 since the fixed point is real and satisfies the real-to-real function needed
c1(x,b)=exp(x*L(b));
c2(x,b)=((exp(x*L(b))*(exp(x*L(b))-1))/(exp(L(b))-1));
\\ here k basically recenters the superfunction so that superf(0)=k.
superfap(x,b,k) = L(b) + c1(x,b)*(k-L(b)) + (c2(x,b)/2!)*((k-L(b))^2);
\\ This recursion is needed to approximate the superfunction of f(x) = exp(x) + b better. its a numerical trick
r(x, b, n) = if (n<-1,log(r(x,b,n+1)-b),if(n==-1,log(x-b),if(n==0,x,if(n==1,exp(x)+b,if(n>1,exp(r(x,b,n-1))+b)))));
\\ the n for superfimp(x,b,n) recenters the superfunction superf, then recenters it back so that superfimp(0,b,n) = 1, this is to improve accuracy
superfimp(x,b,n) = r(superfap(x,b,r(1,b,-n)),b,n);
\\takes a taylor series centered at -1.05 to calculate the superfunction f(x) = exp(x) + b where b=0 so we can just get sexp(x)
\\ This code has an avg error of around ~10^-15 on the specified interval. its pretty good. the only limiting factor here is how we get from b = -1.05 to b = 0
sexp(o) = {
default(realprecision, 300);
my(L = subst(bestapprPade(Ser(superfimp(o,x-1.05,300),x,30)),x,1.05));
return(precision(L,38));
}Let me explain the general process of how I get a superfunction. Also, I would like to mention that the things I mention here may or may not be accurate, I use the things i am explaining here to numerically solve for superfunctions, this is not a rigorous framework, it was built to compute.
Let \( \operatorname{g}(x) \) be the Super function of \( \operatorname{f}(x) \) which means that \( \operatorname{g}(x) = \operatorname{f}(\operatorname{g}(x-1)) \). Let \(\operatorname{g}(0) = k\).
For \(n \in \mathbb{C}\):
\[
\operatorname{g}(x) = \operatorname{f}^{\circ (n)}(\operatorname{g}(x-n))
\]
Using \( \operatorname{g}(x) = \operatorname{f}^{\circ (n)}(\operatorname{g}(x-n)) \) we can then:
\[
\operatorname{g}(x) = \operatorname{f}^{\circ (n)}(\operatorname{g}(x-n)) \Rightarrow \operatorname{g}(x+n) = \operatorname{f}^{\circ (n)}(\operatorname{g}(x))
\]
And so substituting \(x\) with \(\operatorname{g}^{\circ (-1)}(x)\) in \(\operatorname{g}(x+n) = \operatorname{f}^{\circ (n)}(\operatorname{g}(x)) \) we get that:
\[
\operatorname{f}^{\circ (n)}(x) = \operatorname{g}(\operatorname{g}^{\circ (-1)}(x)+n)
\]
We know that \(\operatorname{g}(0) = k\) and so applying the function \(\operatorname{g}^{\circ (-1)}(x)\) on both sides for \(\operatorname{f}^{\circ (n)}(x) = \operatorname{g}(\operatorname{g}^{\circ (-1)}(x)+n)\) we get:
\[
\operatorname{g}^{\circ {-1}}(k) = 0
\]
And so substituting \(x\) with \(k\) in \(\operatorname{f}^{\circ (n)}(x) = \operatorname{g}(\operatorname{g}^{\circ (-1)}(x)+n) \) we get for \(n \in \mathbb{C}\):
\[
\operatorname{f}^{\circ (n)}(k) = \operatorname{g}(n)
\]
This dissolves the problem into finding the nth functional iteration of \(\operatorname{f}(x)\). Now we can start solving for a local solution for \(\operatorname{f}^{\circ (n)}(x)\) around the fixed point of \(\operatorname{f}(x)\).
_________________________________________________________________________________________________________________________________________________________________________________________
Let \(L=\operatorname{f}(L)\) and \(\lambda = \operatorname{f}'(L)\) where \(\lambda\) must not equal \(0\).
For \(n \in \mathbb{C}\) we can actually solve for the first derivative of \(\operatorname{h}(x)\) where \(\operatorname{h}^{\circ (n)}(x) = \operatorname{f}(x)\) which means that \(\operatorname{h}(x) = \operatorname{f}^{\frac{1}{n}}(x)\).
But first we need to establish an identity where \(\operatorname{h}(L) = L\). Now we can continue.
For example for \(n=2\):
\[
\operatorname{h}(\operatorname{h}(x)) = \operatorname{f}(x)
\]
taking the derivative of both sides:
\[
\frac{d}{dx}\operatorname{h}(\operatorname{h}(x)) = \frac{d}{dx}\operatorname{f}(x)
\]
\[
\operatorname{h}'(x) \cdot \operatorname{h}'(\operatorname{h}(x)) = \operatorname{f}'(x)
\]
substituting \(x\) with \(L\):
\[
\operatorname{h}'(L) \cdot \operatorname{h}'(L) = \operatorname{f}'(L)
\]
and so:
\[
\operatorname{h}'(L) = \operatorname{f}'(L)^{\frac{1}{2}} = \lambda^{\frac{1}{2}}
\]
We can do the same for \(n=3\):
\[
\operatorname{h}(\operatorname{h}(\operatorname{h}(x))) = \operatorname{f}(x)
\]
taking the derivative of both sides:
\[
\frac{d}{dx}\operatorname{h}(\operatorname{h}(\operatorname{h}(x))) = \frac{d}{dx}\operatorname{f}(x)
\]
\[
\operatorname{h}'(x) \cdot \operatorname{h}'(\operatorname{h}(x)) \cdot \operatorname{h}'(\operatorname{h}(\operatorname{h}(x))) = \operatorname{f}'(x)
\]
substituting \(x\) with \(L\):
\[
\operatorname{h}'(L) \cdot \operatorname{h}'(L) \cdot \operatorname{h}'(L) = \operatorname{f}'(L)
\]
and so:
\[
\operatorname{h}'(L) = \operatorname{f}'(L)^{\frac{1}{3}} = \lambda^{\frac{1}{3}}
\]
We can do this many times for other \(n\) where \(n \in \mathbb{N}\), and we quickly find a pattern where:
\[
\operatorname{h}'(L) = \lambda^{\frac{1}{n}}
\]
And since \(\operatorname{h}(x) = \operatorname{f}^{\circ (\frac{1}{n})}(x)\) then:
\[
(\operatorname{f}^{\circ (\frac{1}{n})})'(L) = \lambda^{\frac{1}{n}}
\]
and so:
\[
(\operatorname{f}^{\circ (n)})'(L) = \lambda^{n}
\]
To get higher order derivatives we can use a general formula where \(\operatorname{B}_{k,n}\) are the Partial Bell Polynomials and \(\lambda \not= 1\)(i found this formula on some math forum i cant find anymore):
\[
(\operatorname{f}^{\circ (n)})^{c}(L) = \frac{\left(\sum_{m=2}^{c} \operatorname{f}^{(m)}(L) \cdot \operatorname{B}_{c,m}\left( (\operatorname{f}^{\circ (n)})'(L), (\operatorname{f}^{\circ (n)})''(L), \dots, (\operatorname{f}^{\circ (n)})^{(c-m+1)}(L) \right)\right) - \left(\sum_{m=1}^{c-1} \left((\operatorname{f}^{\circ (n)})^{(m)}(L) \cdot \operatorname{B}_{c,m}(\operatorname{f}'(L), \operatorname{f}''(L), \dots, \operatorname{f}^{(c-m+1)}(L))\right) \right)}{\lambda^{c} - \lambda}
\]
This formula seems to work well after manually checking it numerically for different derivatives and functions.
However for \(\lambda = 1\)(parabolic fixed point) the method is quite different since the formula above only works for a fixed point where \(\lambda \not= 1\).
I have made some python code for the fixed point where \(\lambda = 1\):
Code:
# this is a resource, basically for derivative when f'(L)=1. this is Version 56, lets hope i dont need to make more versions.
import math
from fractions import Fraction
def bellpepper(n, k):
# Generates partial Bell polynomial(i forgot how this works again, if it works dont touch it)
if n == 0 and k == 0: return [(1, ())]
if n == 0 or k == 0: return []
res = []
for m in range(1, n - k + 2):
for coeff, indices in bellpepper(n - m, k - 1):
res.append((coeff * math.comb(n - 1, m - 1), (m,) + indices))
return res
def evaluate_sequence(n_max):
# Simulates the true numerical steps of the derivative histat a fixed point where f(L) L and f'(L) = 1, if not then it simulates something atleast
hist = []
# KEEP IT CLEAN. KEEP IT CLEAN. KEEP IT CLEAN. REMEMBER YOU CANT READ THE CODE WHEN DEBUGGING
h0 = {1: {(): Fraction(1)}}
for i in range(2, n_max + 1): h0[i] = {}
hist.append(h0)
for x in range(1, n_max + 2):
prev = hist[-1]
curr = {1: {(): Fraction(1)}}
for n in range(2, n_max + 1):
curr_n = {}
for a in range(1, n + 1):
bell_terms = bellpepper(n, a)
for b_coeff, indices in bell_terms:
term_profile = {(): Fraction(b_coeff)}
for idx in indices:
new_profile = {}
for p1, c1 in term_profile.items():
for p2, c2 in prev[idx].items():
m = max(len(p1), len(p2))
e1 = p1 + (0,) * (m - len(p1))
e2 = p2 + (0,) * (m - len(p2))
f_res = tuple(a_val + b_val for a_val, b_val in zip(e1, e2))
new_profile[f_res] = new_profile.get(f_res, 0) + c1 * c2
term_profile = new_profile
for p, c in term_profile.items():
if a >= 2:
m = max(len(p), a - 1)
e = p + (0,) * (m - len(p))
e_list = list(e)
e_list[a - 2] += 1
final_p = tuple(e_list)
else:
final_p = p
curr_n[final_p] = curr_n.get(final_p, 0) + c
curr[n] = {k: v for k, v in curr_n.items() if v != 0}
hist.append(curr)
return hist
def get_nth_derivative_expansion(n_max):
#Performs Langrange interpolation, or atleast tries to
hist = evaluate_sequence(n_max)
formulas = {}
x_points = list(range(len(hist)))
for n in range(1, n_max + 1):
all_profiles = set()
for h in hist:
for p in h[n].keys(): all_profiles.add(p)
formulas[n] = {}
for p in all_profiles:
y_points = [h[n].get(p, Fraction(0)) for h in hist]
poly = {0: Fraction(0)}
for i, y in enumerate(y_points):
if y == 0: continue
term = {0: Fraction(1)}
denom = 1
for j in range(len(x_points)):
if j == i: continue
denom *= (i - j)
new_term = {}
for pow, coeff in term.items():
new_term[pow + 1] = new_term.get(pow + 1, 0) + coeff
new_term[pow] = new_term.get(pow, 0) - coeff * j
term = new_term
scale = Fraction(y, denom)
for pow, coeff in term.items():
poly[pow] = poly.get(pow, 0) + coeff * scale
formulas[n][p] = {k: v for k, v in poly.items() if v != 0}
return formulas
def print_nth_expansion(target_n):
# does something i forgot, im writing these comments 3 months later
if target_n < 1: return
all_formulas = get_nth_derivative_expansion(target_n)
terms = []
for p, x_poly in all_formulas[target_n].items():
f_parts = []
for idx, power in enumerate(p):
if power == 0: continue
deriv_order = idx + 2
base_str = "f''" if deriv_order == 2 else "f'''" if deriv_order == 3 else f"f^({deriv_order})"
f_parts.append(f"({base_str}(L))^{power}" if power > 1 else f"{base_str}(L)")
f_label = " * ".join(f_parts) if f_parts else "1"
x_parts = []
for pow, coeff in sorted(x_poly.items()):
coeff_str = str(coeff) if coeff.denominator == 1 else f"({coeff.numerator}/{coeff.denominator})"
if pow == 0: x_parts.append(f"{coeff_str}")
elif pow == 1: x_parts.append(f"{coeff_str}*x")
else: x_parts.append(f"{coeff_str}*x^{pow}")
x_str = " + ".join(x_parts)
terms.append(f"[{f_label}] * ({x_str})")
print(f"c(x, {target_n}) = " + " + ".join(terms))
if __name__ == "__main__":
# change n to get derivative. yay
target_derivative_order = 5
print_nth_expansion(target_derivative_order)Now that we can solve for \((\operatorname{f}^{\circ (n)})^{c}(L)\), we can approximate a local continuation for \(\operatorname{g}(x)\) using \((\operatorname{f}^{\circ (n)})^{c}(L)\):
\[
\operatorname{g}(x) \approx\ L + \sum_{n=1}^{\infty} (\operatorname{f}^{\circ (x)})^{n}(L)(k-L)^{n}
\]
or:
\[
\operatorname{g}(x) \approx\ L + \lambda^{x}(k-L) + \frac{\operatorname{f}''(L)(\lambda^{2x} - \lambda^{x})}{\lambda^{2} - \lambda}(k-L)^{2} ...
\]
I have some important conditions when numerically calculating the local approximation of \(\operatorname{g}(x)\) to ensure stability(this is a work in progress and more conditions maybe added later on):
1. If \( \lim_{n \to -\infty} \operatorname{g}(n) \) diverges and \( \lim_{n \to \infty} \operatorname{g}(n) \) converges and is the fixed point of \( \operatorname{f}(x) \) then that should be L, equivalently if \( \lim_{n \to \infty} \operatorname{g}(n) \) diverges and \( \lim_{n \to -\infty} \operatorname{g}(n) \) converges and is the fixed point of \( \operatorname{f}(x) \) then that should be L. If both \( \lim_{n \to -\infty} \operatorname{g}(n) \) and \( \lim_{n \to \infty} \operatorname{g}(n) \) converges and is the fixed point of \(\operatorname{f}(x)\) then either of them could be used.
2. The region where \(\operatorname{g}(x)\) converges to the fixed point \(L\) must be almost linear, if not then it is going to be a pain trying to use the local approximation for \(\operatorname{g}(x)\).
3. If the function \(\operatorname{g}(x)\) is meant to output real numbers for \(x \in \mathbb{R}\) then the fixed point of \(\operatorname{f}(x)\) should be real and \(\lambda \not= 1\), this condition is important because a real-to-real function with a complex fixed point is hard to force that real-to-real condition on.
4. \(\operatorname{g}(x)\) must not contain a lot of exploding singularities in the direction where the function \(\operatorname{f}(x)\) fixed point is, if it does then you cannot use a finite series and you must use another method to bypass those exploding singularities.
Now that we have the numerical conditions laid out, I can show a numerical technique to get tons of precision out of a finite series for \(\operatorname{g}(x)\), where \(\operatorname{g}(x)\) follows the numerical conditions.
Now lets define a function where \(N\) is a finite real number:
\[
\operatorname{g}_{approx}(x,k) = L + \sum_{n=1}^{N} (\operatorname{f}^{\circ (x)})^{n}(k-L)^{n}
\]
I usually use \(N=2\) which is a quadratic series. Let \(dir\) be the direction where the function \(\operatorname{g}(x)\) converges to the fixed point \(L\)(for example: \(dir = -\infty\)).
If we set \(k\) in \(\operatorname{g}_{approx}(x,k)\) closer to the fixed point we see that the function is almost linear there. We can check how good the approximation is by doing \(\\operatorname{f}(\operatorname{g}_{approx}(0,k)) - \operatorname{g}_{approx}(1,k)\), as \(\operatorname{f}(\operatorname{g}_{approx}(0,k)) - \operatorname{g}_{approx}(1,k)\) gets closer to 0 the better the approximation is. And so a improved approximation of \(\operatorname{g}(x)\) is:
\[
\operatorname{f}^{\circ (-dir)}(\operatorname{g}_{approx}(x,\operatorname{f}^{\circ (dir)}(k)))
\]
This is what I call 'the re-centering trick'
Now that I have got the basics laid out, I can now explain how the code works.
First of all we cannot directly calculate \(\operatorname{g}(x) = e^{\operatorname{g}(x-1)}\) using this method, as it is a function where a large portion of its real plane is meant to be real(\(x \in [-2,\infty]\)) but its fixed point is complex(which violates numerical condition 3), and in the direction where the function converges to the fixed point there is a lot of singularities(which violates numerical condition 4).
So instead of looking at \(\operatorname{g}(x) = e^{\operatorname{g}(x-1)}\), instead for \(b < -1\) we look at:
\[
\operatorname{g}(x) = e^{\operatorname{g}(x)} - b
\]
This means that \(\operatorname{f}(x) = e^{x} - b\)
The reason we look at the function \(\operatorname{g}(x) = e^{\operatorname{g}(x)} - b\) for \(b<-1\) is because the fixed point is real and the function is meant to output a real number for all reals which means it satisfies numerical condition 3, also for \(b < -1\) the function has no singularities in its real plane which satisfies numerical condition 4, and the region where the function \(\operatorname{g}(x)\) converges to the fixed point is almost linear which satisfies numerical condition 2.
According to Numerical condition 1 the fixed point would be:
\[
L = b - \operatorname{W}_{-1}(-e^{b})
\]
Where \(\operatorname{W}_{-1}\) is the -1th branch of the Lambert W function. The reason we are using the -1th branch is because that the direction where \(\operatorname{g}(x)\) converges to the fixed point is to the left of the real axis(meaning that \(dir = -\infty\)) and the -1th branch gets us that fixed point.
In my code the fixed point for \(b\) is defined as \(L(x) = x - \operatorname{W}_{-1}(-e^{x})\) with this i construct a quadratic approximation of \(\operatorname{g}(x) = e^{\operatorname{g}(x)} - b\):
\[
\operatorname{superfap}(x,b,k) = L(b) + e^{x \cdot L(b)}(k-L(b)) + \frac{e^{x \cdot L(b)} \cdot (e^{x \cdot L(b)} - 1)}{e^{L(b)} - 1}
\]
This is not accurate enough at \(k = 1\) so what we do is use our recentering trick from earlier to improve accuracy by a ton:
\[
\operatorname{superfimp}(x,b,n) = \operatorname{f}^{\circ (n)}(\operatorname{superfap}(x,b,\operatorname{f}^{\circ (-n)}(1)))
\]
In my code \(\operatorname{r}(x,b,n) = \operatorname{f}^{\circ (n)}(x)\). To get \(\operatorname{superfimp}(x,0,n)\)(this is just \(\operatorname{sexp}(x)\)) where in my code \(n\) is set to 300 we have to analytically continue from a point close to -1 to 0. In my code the point is hardcoded as -1.05 because at -1 \(\lambda = 1\) which makes the coefficient of \((k-L)^{2}\) equal to 0, and so the approximation becomes unstable near -1, and that is why i use -1.05 instead of something like -1.02. We can just take a taylor series at \(x=-1.05\), which is exactly what V1 does which achieves an accuracy of ~\(10^{-10}\)(V1 uses a 18 degree taylor series because that is the most accurate as it can get):
Code:
\\ Numerical algorithm to calculate an approximation of Kneser tetration on the interval [0,1]
\\ this is V1 which includes taking a series using Ser(expr,var,degree) to make the error near ~10^-10 on the specificied interval
\\ fixed point for f(x) = exp(x) + b
default(seriesprecision, 30);
L(x) = x-lambertw(-exp(x),-1);
\\ faa-di-bruno reccurence series for the superfunction of f(x) = exp(x) + b.
\\ Works for b < -1 since the fixed point is real and satisfies the real-to-real function needed
c1(x,b)=exp(x*L(b));
c2(x,b)=((exp(x*L(b))*(exp(x*L(b))-1))/(exp(L(b))-1));
\\ here k basically recenters the superfunction so that superf(0)=k.
superfap(x,b,k) = L(b) + c1(x,b)*(k-L(b)) + (c2(x,b)/2!)*((k-L(b))^2);
\\ This recursion is needed to approximate the superfunction of f(x) = exp(x) + b better. its a numerical trick
r(x, b, n) = if (n<-1,log(r(x,b,n+1)-b),if(n==-1,log(x-b),if(n==0,x,if(n==1,exp(x)+b,if(n>1,exp(r(x,b,n-1))+b)))));
\\ the n for superfimp(x,b,n) recenters the superfunction superf, then recenters it back so that superfimp(0,b,n) = 1, this is to improve accuracy
superfimp(x,b,n) = r(superfap(x,b,r(1,b,-n)),b,n);
\\takes a taylor series centered at -1.05 to calculate the superfunction f(x) = exp(x) + b where b=0 so we can just get sexp(x)
\\ This code has an avg error of around ~10^-10 on the specified interval. its pretty good. the only limiting factor here is how we get from b = -1.05 to b = 0
sexp(o) = {
default(realprecision, 300);
my(L = subst(Pol(Ser(superfimp(o,x-1.05,300),x,18)),x,1.05));
return(precision(L,38));
}However in V2 it uses a 30 degree taylor series then applies a Pade approximant to it to accelerate convergence making the precision higher(around ~\(10^{-15}\)). Also to increase precision, you would need to use a better numerical technique to analytically continue from \(b=-1.05\) to \(b=0\). I know I didn't provide numerical error analysis in this post, that is because im kinda lazy to do it. Maybe in the near future I could provide some numerical error analysis for V2. Anyways for V3 i might generalize the algorithm for all bases except for 0 and 1 and maybe in a later version use a Borel-conformal resumption. This was a fun side project to do. If you do have any questions, feel free to ask!