Tetration Forum - All Forums

PARI/GP code to approximate Kneser tetration(base e) without using schroder/abel

Mon, 27 Jul 2026 16:45:21 +0000

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)

And the code does seem to work very well. I have used this code to extend Product Tetration: https://tetrationforum.org/showthread.php?tid=1808.
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} \frac{(\operatorname{f}^{\circ (x)})^{n}(L)}{n!}(k-L)^{n}
\]

or:

\[
\operatorname{g}(x) \approx\ L + \lambda^{x}(k-L) + \frac{\operatorname{f}''(L)(\lambda^{2x} - \lambda^{x})}{2(\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)}{2(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!

Mixed-base tetration: five preprints, a faster fatou.gp fork, and a base-change calcu

Wed, 15 Jul 2026 07:33:05 +0000

Hello everyone,
my name is Janis, and I am a high-school student with a strong interest in tetration. Over the last one or two years, I have spent much of my free time experimenting with tetration numerically and trying to understand some of its underlying structure.
The work did not begin as a rigorous research project. It started with computations, observations, and many incomplete ideas. With the help of AI tools, I was able to test these ideas, formulate them more precisely, and develop them into a series of papers. The papers contain a detailed statement explaining the role of AI in the project, and I take responsibility for the final content.
Today I am sharing three connected parts of the project:
1. A five-paper series on mixed-base tetration
The main object I study is the mixed-base height map
F_{b,d}(x) = slog_b(T_d(x)),
where T_d is tetration in base d and slog_b is the superlogarithm in base b.
This provides a type of base change between different tetration functions. Numerically, the difference between this map and the original height does not generally approach a constant. Instead, it appears to approach a small periodic phase function. The papers investigate this phase, its algebraic structure, its interpretation as a circle map, its numerical fine structure, and its dependence on the tetration base.
All five papers are available as preprints on Zenodo, and their PDFs are included in the repository:
https://github.com/Lightrunnerwastaken/g...ain/papers
2. An optimized fork of Sheldonison’s fatou.gp
All computations build on Sheldonison’s original fatou.gp, for which I am very grateful.
I made a number of performance improvements while repeatedly comparing the results against frozen accuracy tests. On my machine, one high-precision base-e computation at approximately 500 digits was reduced from about 5.4 hours to about 25 minutes. The exact benchmarks and methods are documented in the repository.
The repository also contains high-precision reference values with numerical error bounds. The strongest current result is a computation of sexp_e(0.5) with a claimed 972-digit error bound.
3. A base-change calculator
Using the phase functions, I constructed a small “base atlas.” After preparing one anchor tetration function, it approximates sexp and slog for other bases using a short table of Fourier coefficients..
The complete repository is here:
https://github.com/Lightrunnerwastaken/gp-tetration
I would be very grateful for feedback, especially about possible errors, missing assumptions, the claimed error bounds, and the mathematical interpretation of the phase functions. I am still learning and would be happy to correct or improve anything.
I also have a practical question for Sheldonison: I could not find an explicit license in fatou.gp. I have preserved the attribution and documented the origin of the code. Please let me know whether you are comfortable with the fork being distributed in this form and what licensing terms you would prefer. If you would rather have the original file removed from the repository, I will of course do so.

Preprints:

A Phase Law for Mixed-Base Tetration
https://doi.org/10.5281/zenodo.21346803
The Mixed-Base Tetration Phase Map Is a Circle Diffeomorphism
https://doi.org/10.5281/zenodo.21361967
Stopped Channels and Offset Alignment
https://doi.org/10.5281/zenodo.21362206
Fine Structure of Mixed-Base Tetration Phases
https://doi.org/10.5281/zenodo.21363593
The Base-Derivative of the Kneser Family
https://doi.org/10.5281/zenodo.21363809

Seeking an arXiv endorser for mixed-base tetration (math.DS)

Mon, 13 Jul 2026 21:19:29 +0000

Hello,
I am a student and independent researcher working on a broader research project concerning mixed-base tetration. I am currently preparing the first manuscript from this project for submission to arXiv:
“A Phase Law for Mixed-Base Tetration: Logarithmic Tails, Cocycle Structure, and Singular Asymptotics.”
The paper studies mixed-base tetration through Abel coordinates, logarithmic tails, periodic phase maps, cocycle identities, and asymptotic behavior. Several related manuscripts are also in preparation and are intended to be released sequentially after the first paper.
As a first-time arXiv submitter, I require an endorsement for math.DS. Is there an established arXiv author here who would be willing to assess whether the manuscript is scholarly and appropriate for this category?
I would be happy to provide the complete PDF privately.

My introduction and Some Approximations of Kneser Tetration on Desmos I made

Tue, 07 Jul 2026 13:22:44 +0000

Hello! My name is Raees. I have been researching tetration and hyperoperations for 5 months. I am currently studying in High School. My birthday is July 2nd. I was born in India. My Hobby is currently researching Tetration and Hyperoperations(since there is so much to explore!) and making numerical models of functions in Desmos. I am mostly self-taught in math.

Anyways, I have tried to make some implementations of Kneser Tetration for Complex heights and a fixed base on Desmos(specifically the function sexp(a+bi), with the data computed using PARI/GP and the program fatou.gp made by sheldonism):

https://www.desmos.com/calculator/hsjtrcgqgs (Test Version to check my idea, works for imag(x) ∈ [-5,5], base e)

These are the Current Versions for a variety of Bases(which technically works for all complex numbers, but i reccommend the real part staying low ):
https://www.desmos.com/calculator/n1prudguef (Current version, base 1.5)
https://www.desmos.com/calculator/anektsk1pb (Current version, base 2)
https://www.desmos.com/calculator/d4frqmds2e (Current version, base 2.5)
https://www.desmos.com/calculator/otrjcsnhoe (Current version, base e)
https://www.desmos.com/calculator/rwfsrv1lte (Current version, base 3)
https://www.desmos.com/calculator/dtvyzocshm (Current version, base Pi)
https://www.desmos.com/calculator/8jbdph3pgv (Current version, base i(I'm not actually sure if its Kneser, i got the data from fatou.gp made by sheldonison))

I also made another approximation for real heights and bases ∈ [1.5,3.5]: https://www.desmos.com/calculator/12qaelj9vz(data computed from fatou.gp)

_________________________________________________________________________________________________________________________________________________

My current goal is to make Kneser Tetration for a lot of bases in Desmos. I discovered a path were I think i could achieve an approximation of Kneser tetration(but I need to research and refine my method, also a little sneak peak into the method: https://www.desmos.com/calculator/wonk57usxr).
Feel free to ask me any questions about my approximations of Kneser Tetration on Desmos!(also I am known as DubiumTetration on Reddit, and i uploaded one of my graphs on reddit:"https://www.reddit.com/r/desmos/comments/1tzgjnl/kneser_tetration_for_base_e_for_complex_heights/")

test

Wed, 24 Jun 2026 17:38:14 +0000

\[a+b\], \(c+d\)

[tikz2]\draw (0,0) circle (1in);[/tikz2]

[tikz2]A\arrow[r, "a"]& B\\ [/tikz2]

can normal tetration be defined as a limit of product tetration?

Wed, 27 May 2026 17:00:54 +0000

it is hypothesized that if you take infinitely many arbitrary base logarithms of product tetration with the same base, you will eventually get normal tetration with the base you chose.
but is it really true? someone please let me know.
here's an image to make things clearer:

image of tetration.png (Size: 85.42 KB / Downloads: 215)

Recurrence relations and differential equations

Wed, 06 May 2026 06:58:49 +0000

I'd like to hear your thoughts on why recurrence relations are considered the discrete version of differential equations.
I'm trying to generalize recurrence relations and autonomous differential equations.
I think that \(\mathbb{R}_+\)-actions are the straightforward continuous version of \(\mathbb{N}\)-actions, so I checked whether an orbit of a \(\mathbb{R}_+\)-action is always a solution (an integral curve) of an autonomous differential equation, and the opposite.

\(\gamma: \mathbb{R} \to \mathbb{R}^n\) is an orbit of a \(\mathbb{R}\)-action if and only if \(\gamma(s) = \gamma(t) \implies \gamma(s+u) = \gamma(t+u)\) for all \(u \in \mathbb{R}\).
\(\gamma\) is an orbit of a \(\mathbb{R}_+\)-action if and only if \(\gamma(s) = \gamma(t) \implies \gamma(s+u) = \gamma(t+u)\) for all \(u > 0\).
\(\gamma\) is a solution of an autonomous differential equation \(\gamma'(t) = F(\gamma(t))\) if and only if \(\gamma(s) = \gamma(t) \implies \gamma'(s) = \gamma'(t)\). (Define \(F(x) := \gamma'(t)\) when \(x = \gamma(t)\). The \(F\) is well-defined on \(\gamma(\mathbb{R})\) because of the hypothesis.)

From the Picard–Lindelöf theorem, if \(F: \mathbb{R}^n \to \mathbb{R}^n\) is Lipschitz continuous, a solution \(\gamma: \mathbb{R} \to \mathbb{R}^n\) of a differential equation \(\gamma'(t) = F(\gamma(t))\) is always an orbit of a \(\mathbb{R}\)-action.
Trivially an orbit of a \(\mathbb{R}\)-action is an orbit of a \(\mathbb{R}_+\)-action.
If an orbit of a \(\mathbb{R}_+\)-action is differentiable, it's always a solution of an autonomous differential equation. (The proof is described below.)
Therefore let \(\gamma: \mathbb{R} \to \mathbb{R}^n\) is differentiable, then
\[
\begin{array}{ll}
& \text{there exists a Lipschitz continuous \(F: \mathbb{R}^n \to \mathbb{R}^n\) such that \(\gamma'(t) = F(\gamma(t))\)} \\
\implies & \text{\(\gamma\) is an orbit of a \(\mathbb{R}\)-action} \\
\implies & \text{\(\gamma\) is an orbit of a \(\mathbb{R}_+\)-action} \\
\implies & \text{there exists \(F: \mathbb{R}^n \to \mathbb{R}^n\) such that \(\gamma'(t) = F(\gamma(t))\)}.
\end{array}
\]
The opposites are not true.

An orbit of a \(\mathbb{R}\)-action that is not a solution of an autonomous differential equation of a Lipschitz continuous function:
Since
\[
\gamma(t) =
\begin{cases}
-t^2 & (t \leq 0) \\
t^2 & (t \geq 0)
\end{cases}
\]
is injective, the \(\gamma: \mathbb{R} \to \mathbb{R}^1\) is an orbit of \(\mathbb{R}\)-action.
Assume there exists a Lipschitz continuous \(F: \mathbb{R}^1 \to \mathbb{R}^1\) such that \(\gamma'(t) = F(\gamma(t))\).
Then there exists \(L \in \mathbb{R}\) such that \(|\gamma'(t_1) - \gamma'(t_2)| = |F(\gamma(t_1)) - F(\gamma(t_2))| \leq L|\gamma(t_1) - \gamma(t_2)|\) for all \(t_1, t_2 \in \mathbb{R}\).
Let \(t_1 \neq 0\) and \(t_2 = 0\).
We see \(\gamma'(t) = 2|t|\), so \(2|t_1| \leq L|t_1^2|\). Therefore
\[
\left|\frac{1}{t_1}\right| = \left|\frac{t_1}{t_1^2}\right| \leq \frac{L}{2}.
\]
\(\left|\frac{1}{t_1}\right| \longrightarrow \infty\) when \(t_1 \longrightarrow 0\), that is a contradiction.

A differentiable orbit of a \(\mathbb{R}_+\)-action that is not an orbit of a \(\mathbb{R}\)-action:
\[
\gamma: \mathbb{R} \to \mathbb{R}^1, \quad
\gamma(t) =
\begin{cases}
t^2 & (t \leq 0) \\
0 & (t \geq 0)
\end{cases}
\]
is an example of it.

A solution of an autonomous differential equation that is not an orbit of a \(\mathbb{R}_+\)-action:
\[
\gamma: \mathbb{R} \to \mathbb{R}^1, \quad
\gamma(t) =
\begin{cases}
0 & (t \leq 0) \\
t^2 & (t \geq 0)
\end{cases}
\]
is a solution of \(\gamma'(t) = F(\gamma(t))\) where
\[
F: \mathbb{R}^1 \to \mathbb{R}^1, \quad
F(x) =
\begin{cases}
0 & (x \leq 0) \\
2\sqrt{x} & (x \geq 0),
\end{cases}
\]
while \(\gamma(-1) = \gamma(0)\) and \(\gamma(-1+1) = 0 \neq 1 = \gamma(0+1)\).

Proof of that a differentiable orbit of a \(\mathbb{R}_+\)-action is a solution of an autonomous differential equation:
We'll prove the contraposition.
Let \(\gamma: \mathbb{R} \to \mathbb{R}^n\) is differentiable, \(\gamma(t) = (\gamma_1(t), \dots, \gamma_n(t))\). (\(\gamma'(t)\) denotes \((\gamma_1'(t), \dots, \gamma_n'(t))\).)
Assume \(\gamma\) is not a a solution of an autonomous differential equation.
Then there are \(t_1, t_2 \in \mathbb{R}\) such that \(\gamma(t_1) = \gamma(t_2)\) and \(\gamma'(t_1) \neq \gamma'(t_2)\).
Therefore \(\gamma'_i(t_1) \neq \gamma'_i(t_2)\) for some \(i\).
Define \(x = \gamma_i(t_1) = \gamma_i(t_2)\), \(\alpha = \gamma_i'(t_1)\), and \(\beta = \gamma_i'(t_2)\). We can assume \(\alpha < \beta\) w.l.o.g.
We know
\[
\forall \varepsilon_1 > 0, \exists \delta_1 > 0, \forall h \in (-\delta_1, \delta_1), \left|\frac{\gamma_i(t_1+h)-x}{h} - \alpha\right| < \varepsilon_1
\]
and
\[
\forall \varepsilon_2 > 0, \exists \delta_2 > 0, \forall h \in (-\delta_2, \delta_2), \left|\frac{\gamma_i(t_2+h)-x}{h} - \beta\right| < \varepsilon_2.
\]
Thus
\[
x + \alpha h - \varepsilon_1|h| < \gamma_i(t_1+h) < x + \alpha h + \varepsilon_1|h|,
\]
\[
x + \beta h - \varepsilon_2|h| < \gamma_i(t_2+h) < x + \beta h + \varepsilon_2|h|.
\]
Since \(\beta - \alpha > 0\) there are \(\varepsilon, \varepsilon' \in \mathbb{R}\) such that \(0 < \varepsilon < \varepsilon' < \beta - \alpha\).
Substitute \(\varepsilon_1 = \varepsilon\) and \(\varepsilon_2 = \beta - \alpha - \varepsilon'\), then for any \(0 < h < \min\{\delta_1, \delta_2\}\),
\begin{align*}
\gamma_i(t_1+h) &< x + \alpha h + \varepsilon_1|h| \\
&= x + \alpha h + \varepsilon h \\
&< x + \alpha h + \varepsilon'h \\
&= x + \beta h - (\beta - \alpha - \varepsilon')h \\
&= x + \beta h - \varepsilon_2|h| \\
&< \gamma_i(t_2+h).
\end{align*}
That is \(\gamma(t_1+h) \neq \gamma(t_2+h)\).
Therefore \(\gamma\) is not an orbit of a \(\mathbb{R}_+\)-action.

From these results, I lost sight of what is the relation between an autonomous differential equation and a \(\mathbb{R}_+\)-action (or a \(\mathbb{N}\)-action).
I was hoping you could help me find that relation.

the fraction iteration approximation of tetration

Tue, 14 Apr 2026 17:52:40 +0000

we approximate the arbitrary base logarithm with a rational function, namely (ax+b)/(cx+d).
we choose a to be 2, b to be -2, c to be ln(a1), d to also be ln(a1)
a1 is the arbitrary base of the logarithm.
when we're done, we use the iteration formula for (ax+b)/(cx+d).
and we get this graph: https://www.desmos.com/calculator/x6malmxjzt

Reviewing recent developments?

Wed, 31 Dec 2025 06:09:33 +0000

It’s been a long time since I participated in community discussions. Have there been any significant new developments regarding tetration lately? Also, has the community discussed the paper Holomorphic Extension of Tetration to Complex Bases and Heights via Schröder's Equation (https://doi.org/10.13140/RG.2.2.10348.48008)? What are its primary conclusions?
How this method compares to the Kneser?

extending normal tetration to real numbers using product tetration

Fri, 12 Dec 2025 18:49:13 +0000

we start by constructing the base e product tetration of x which is denoted as \(\text{prodtet}_{e}(x)\)
the base e product tetration has these recursive properties:
\(\text{prodtet}_{e}(x)=\text{prodtet}_{e}(x-1)*e^{\text{prodtet}_{e}(x-1)}\)
\(\text{prodtet}_{e}(x)=\text{W}(\text{prodtet}_{e}(x+1))\)
also, product tetration has this following graph: https://www.desmos.com/calculator/zkjiz8jwh0
if you're wondering how i extended product tetration to real numbers, i used lagrange interpolation on the very far negative side of product tetration and applied recursive formulas.
some example values of product tetration:
\(\text{prodtet}_{e}(-1)=\text{W}(1)=0.56714... \)
\(\text{prodtet}_{e}(0)=1\)
\(\text{prodtet}_{e}(1)=e\)
\(\text{prodtet}_{e}(2)=e^{e+1}\)
and one accepted non-integer value for product tetration:
\(\text{prodtet}_{e}(1/2)=1.5134...\)
and we can express product tetration to different bases by introducing the inverse function of product tetration which will be called the product super logarithm.
the product super logarithm will be denoted as \(\text{prodslog}_{e}(x)\)
the product super logarithm also has some noteworthy recursive properties of its own:
\(\text{prodslog}_{e}(x)=\text{prodslog}_{e}(xe^{x})-1\)
\(\text{prodslog}_{e}(x)=\text{prodslog}_{e}(\text{W}(x))+1\)
the product super logarithm also has a noteworthy graph (sorry if it's slow): https://www.desmos.com/calculator/sxjnv2212j
now for the base change formula of both the product tetration and the product super logarithm:
\( \text{prodtet}_{a}(x)=\frac{\text{prodtet}_{e}(\text{prodslog}_{e}(\text{ln}(a))+x)}{\text{ln}(a)}\)
\(\text{prodslog}_{a}(x)=\text{prodslog}_{e}(x*\text{ln}(a))-\text{prodslog}_{e}(\text{ln}(a))\)
here, i give an example of the base change formulas with base pi: https://www.desmos.com/calculator/kwsqnpe4ef
the proof of the base change formulas is left as an exercise to the reader.
now, back to the main topic of this post.
we can literally extend normal tetration to real numbers by using product tetration because product tetration grows faster than normal tetration.
here's how the process goes:
we extend both the product tetration and the product super logarithm to real numbers using any interpolation technique we consider good. (in this case, i considered the lagrange interpolation good.)
we take the product super logarithm of normal tetration and interpolate the result we get by doing so. (the interpolation still has to be good.)
we finally take the product tetration of the above interpolation to get normal tetration.
here's an example of this process being applied to base 2 tetration: https://www.desmos.com/calculator/r26svg3olm
the numerical accuracy on desmos is bad for product tetration and the product super logarithm.
and desmos can only calculate numbers up to \(2^{1024}\).
i hope there is a calculator that functions like desmos but with greater numerical accuracy and the ability to calculate greater numbers than \(2^{1024}\).

Confirmed: the constancy of the congruence speed holds in each squarefree radix

Mon, 01 Dec 2025 19:28:32 +0000

It has been a long time since my last post here, so I hope you will forgive the sudden appearance after one year of silence.
I would like to share a very compact one-page note that summarizes six fundamental equations describing the constant congruence speed of integer tetration bases (greater than 1 and not a multiple of the selected radix) in arbitrary numeral systems \(r > 1\).
These few equations condense the key results established in our referenced recent works, including the complete formulas for prime radices, the explicit radix-6 and radix-10 cases, and a couple of elegant identities valid for all squarefree \(r > 2\).

In particular, I am especially fond of identity (4). It highlights the structural role of the \(r\)-adic fixed point \(1_r\): it considers 5 different parameters on the left-hand side, including the height of the power tower, which only needs to be greater than 1, and for \(k = 0\) (by Mihăilescu's theorem) it returns only perfect powers of exact order \(c\) whose radix-\(r\) constant congruence speed is exactly \(t\). Furthermore, it includes hyper-1 as sum/difference, hyper-2 as product, hyper-3 as exponentiation, and finally the function \(V\) arises from tetration (i.e., hyper-4).

I am attaching these formulas as a figure; the one-page note is also available on Zenodo: Six fundamental equations for the constant congruence speed in any numeral system.

Comments, questions, and suggestions are of course welcome... it may also be interesting to discuss possible extensions or implications for related r-adic fixed-point structures or iterated-exponentiation phenomena (I have also developed a very rough key-exchange idea based on solving fundamental equations of the form \(y^t = y\) for sufficiently large odd values of \(t := t\)(\(r\)) in the ring of \(r\)-adic integers \(\mathbb{Z}_r\)).

Six fundamental equations (reduced).jpg (Size: 228.37 KB / Downloads: 383)

what we should call the hyperoperations in different languages

Sat, 15 Nov 2025 15:12:15 +0000

we already have the english names for hyperoperations: tetration, pentation, hexation, etc. but here are a few suggestions for what we should call the hyperoperations in different languages.
romanian: tetrație, pentație, hexație, etc.
icelandic: 4 rekstur, 5 rekstur, 6 rekstur, etc.
note: in icelandic, rekstur means operation.
german, french: same as english.
second note: we already have the russian, turkish, spanish, and italian names for hyperoperations.
mongolian: 4 үйл, 5 үйл, 6 үйл, etc.
third note: in mongolian, үйл means operation.
hungarian: 4 művelet, 5 művelet, 6 művelet, etc.
fourth note: in hungarian, művelet means operation.
you can add to the list if you know a language other than the ones mentioned here.

Tetration with complex bases

Thu, 13 Nov 2025 09:33:54 +0000

Seeing as how, for real bases, tetration is usually defined by computing the inverse Abel function (α^-1(x)), and mapping the region above the main line where Im(α^-1(x)) = 0 (also extended in the negative direction, where it would otherwise be complex) to the upper-half plane, I was wondering how this would be done for complex bases. Is complex-based tetration computed by just extrapolating from the coefficients of real tetration, or just by using the inverse Abel function, or a different method entirely?

Orbit-like maps on linearly ordered groups

Sat, 08 Nov 2025 07:47:16 +0000

I want to extend the indexes of \(a\) of a recurrence relation \(a_{n+1} = f(a_n)\) to arbitrary linearly ordered groups.
For a map \(a: \mathbb{Z} \to X, n \mapsto a_n\), there exists \(f: X \to X\) such that \(a_{n+1} = f(a_n)\) if and only if \(a_m = a_n \implies a_{m+k} = a_{n+k}\) for any positive integer \(k\).
When \(a: \mathbb{Z} \to X\) satisfies the above condition, the periods of \(a\), \(P(a) = \{m-n \mid a_m = a_n\}\) is a subgroup of \(\mathbb{Z}\).

To generalize to arbitrary linearly ordered groups, I will introduce some definitions.
For a map \(a\) defined on a linearly ordered group \(G\),
\(a\) is peep if \(a(g) = a(h) \implies a(fg) = a(fh)\) for all \(f > 1\).
\(a\) is bib if \(a(g) = a(h) \implies a(gf) = a(hf)\) for all \(f > 1\).
\(a\) is peep bib if \(a\) is peep and bib.
I define periods of \(a\) as
\begin{align}
P(a) &= \{g^{-1}h \mid a(g) = a(h)\} \\
&= \{p \in G \mid \exists g, h \in G, p = g^{-1}h\: \text{and}\: a(g) = a(h)\} \\
&= \{p \in G \mid \exists g \in G, a(g) = a(gp)\}.
\end{align}
A subgroup \(H\) of \(G\) is said to be fee if \(g^{-1}Hg \subset H\) for all \(g > 1\).

Let \(G\) be a bi-ordered group and \(a\) be a map defined on \(G\). (Edit: \(G\) need not to be bi-ordered, just need to be linearly ordered.) Then
\begin{equation}
\begin{array}{ccccc}
& & a\: \text{is}\: \textit{peep bib} & \implies & a\: \text{is}\: \textit{peep} \\
& & \Downarrow & & \Downarrow \\
P(a)\: \text{is normal} & \implies & P(a)\: \text{is}\: \textit{fee} & \implies & P(a)\: \text{is a subgroup} \\
\end{array}
\end{equation}
and the converses of these are not true.

For example, let \(F\) be the free group generated by \(a, b\) and bi-ordred in the way of Definition 3 of https://doi.org/10.4153/CJM-2003-034-2
Then the subgroup \(\langle a \rangle\) is not fee.
The fee subgroup \(\langle b^{F_+} \rangle\) is not normal (shown in the attached pdf
momo (4).pdf (Size: 223.12 KB / Downloads: 272) ).
We define a map \(a\) on \(\mathbb{Z}\) as
\begin{equation}
\dots, -7, -6, -5, 22, -3, 22, -1, 0, 1, 0, 1, 0, 1, 0, 1, \dots
\end{equation}
Then \(a\) is neither peep nor bib, and \(P(a) = 2\mathbb{Z}\).

Is my tetration-webspace still locked?

Mon, 20 Oct 2025 15:28:21 +0000

Dear friends -

it is about 6 weeks, that some correspondent from USA told me he cannot access the files on my website. It seemed so as well from Scotland - but the access from Germany seems to work (I can access everything by ftp:// and even by https:// ).
Now the admin of our small network https://helms-net.de tells me he had tried to reverse some changes he had made in spring, and I would like to check, whether international folks can now access my work again.

Could someone please try to access my math-homepage

        https://go.helms-net.de/math/tetdocs/index.htm

and if that is possible, then to download just a single *.pdf-file as test?

Please tell me the result here, or even better, mail me at

     mailto:gottfried.helms@t-online.de

(possibly include the error-message as well)

Thank you in advance -

Gottfried.