Implement Level-Index Addition

Question

Level-index is a system of representing numbers which has been proposed as an alternative to floating-point formats. It claims to virtually eliminate overflow (and underflow, in its symmetric form) from the vast majority of computations.

The usual floating-point format internally represents a number \$x\$ in terms of a mantissa \$m\in[1,2)\$ and an integer exponent \$e\$, where \$x=m\times2^e\$. In contrast, level-index represents a number by an integer level \$l\ge0\$ and an index \$i\in[0,1)\$, defined by
\$l+i=\psi(x)=\begin{cases} x,&x<1\\ 1+\psi(\log x),&\text{otherwise.} \end{cases}\$
Equivalently, \$l\$ and \$i\$ together represent the power tower number \$x=\psi^{-1}(l+i)=\operatorname{exp}^l(i)\$. Accordingly, \$\psi^{-1}\$ grows very fast; \$\psi^{-1}(5)\$ has more than a million digits to the left of the decimal point.

Task

Given two numbers \$\psi(a),\psi(b)\$, compute their level-index sum \$\psi(a+b)\$. You may choose to input or output the level and index together or separately.

Your solution should be accurate to at least six digits after the decimal point and terminate in a reasonable time for all \$\psi(a+b)<8\$, and may not assume your floating-point representation has infinite precision.

Test cases

ψ(a)    ψ(b)    ψ(a+b)              a           b           a+b
0       0       0.0000000000000000  0           0           0
1       1       1.6931471805599453  1           1           2
2       2       2.5265890341390445  2.718281828 2.718281828 5.436563657
3       3       3.2046791426805201  1.5154262e1 1.5154262e1 3.0385245e1
4       4       4.0161875057657443  3.8142791e6 3.8142791e6 7.6285582e6
5       5       5.0000000044114734  2.3e1656520 2.3e1656520 4.7e1656520
1.70879 2.44303 2.6490358380192970  2.031531618 4.746554831 6.778086449
0.8001  3.12692 3.1367187876590184  0.8001      2.2470167e1 2.3270267e1
2.71828 3.14159 3.2134989633416356  7.774915755 2.368508064 3.1459996e1
4.6322  4.6322  4.6322790036002559  1.79455e308 1.79455e308 3.58910e308
4.80382 4.80382 4.8038229252940133  5.9037e4932 5.9037e4932 1.1807e4933

Answer 1

SageMath, 93 bytes

lambda x,y:h(p(x)+p(y))
h=lambda x:x>1 and 1+h(log(x))or x
p=lambda x:x>1 and exp(p(x-1))or x

Try it online!

Answer 2

Python NumPy, 178 bytes

def f(I):
 i,j=sort(I)
 if 4+i>j<5.5:
  j=[add,logaddexp][l:=i>1](p(j-l),p(i-l))
  while j>=1:l+=1;j=log(j)
  j+=l
 return j
p=lambda x:x>1 and exp(p(x-1))or x
from numpy import*

Attempt This Online!

Saved a bit by borrowing @Noodle9's p function.

Previous Python NumPy, 202 bytes

def f(I):
 i,j=sort(I);l,*c=0,
 if 4+i>j<5.5:
  l=i>1
  for y in i,j:
   r=y%1
   while y>=1+l:r=exp(r);y-=1
   c+=r,
  j=[add,logaddexp][l](*c)
  while j>=1:l+=1;j=log(j)
 return l+j
from numpy import*

Attempt This Online!

I'm sure this can be golfed quite a bit.

How

With the given accuracy requirements we need something to bridge the gap between small values (<4.5, say) where we can use exp without overflowing and large (>5.5) where we can simply use the max. Luckily, NumPy has logaddexp built in which saves us one level of exponentiating.