Calculating 3^3^3^3 (very large exponent / how did Wolfram do it?)

Question

I can't find a right algorithm / struct to calculate the number simply in C or Go. However, a class can easily be created in Python.

At first glance, the calculation seems to be very straight forward. However, when you look at the sample calculation from Wolfram Alpha.

https://www.wolframalpha.com/input/?i=3%5E3%5E3%5E3

This breaks both long long (integer, 18-19 digits) and double (float / IEEE 754, up to e+308 digits, with 17 digits' precision).

However, I can cheat a little with Python, as it will automatically allocate more bytes for integer.

Still, 3^(7.625e+13) takes abnormally very long time... (3^3^3 = 7.625e+13).

import math
from decimal import Decimal


class Int:
    _first = ""
    _last = ""
    _length = None  # Int

    val: int = None  # actual int, if applicable

    def __init__(self, val: int = 0) -> None:
        if isinstance(val, Int):
            if val.val is None:
                self._first = val._first
                self._last = val._last
                self._length = val._length
                return
            self.val = val.val
        else:
            self.val = val

        try:
            float(self.val)
        except OverflowError:
            self._first = self.first
            self._last = self.last
            self._length = self.length

            self.val = None

    @property
    def first(self) -> str:
        if self._first:
            return self._first

        return str(self.val)[:8]

    @property
    def last(self) -> str:
        if self._last:
            return self._last

        return str(self.val)[-8:]

    @property
    def length(self):
        if self._length:
            return self._length

        return Int(len(str(self.val)))

    def exp3(self):
        return Int(3) ** self.val

    def tetrate3(self, n: int):
        first = Int(self)
        for _ in range(n - 1):
            first = first.exp3()

        return first

    def __repr__(self) -> str:
        if self.val is None:
            return f"{self.first}...{self.last} ({self.first[0]}.{self.first[1:]}e+{self.length})"

        return f"{self.val}"

    def __pow__(self, _other):
        base = Int(self)
        exp = Int(_other)

        if base.val and exp.val:
            try:
                float(base.val) ** exp.val
                return Int(base.val ** exp.val)
            except OverflowError:
                pass

        log = Decimal(exp.val) * Decimal(math.log10(base.val))
        fl = math.floor(float(log))

        out = Int()
        out._first = f"{(10 ** float(log - fl)):.7f}".replace(".", "")
        out._last = str(pow(int(base.last), exp.val, 10_000_000_000))[-8:]
        out._length = Int(fl)
        out.val = None

        return out


if __name__ == "__main__":
    # After the third digits may be imprecise
    # => 12579723...00739387 (1.2579723e+3638334640024)
    print(Int(3).tetrate3(4))

Answer 1

Wolfram Alpha is giving you an approximate answer, which is much easier than calculating an exact answer. Most likely it's using the transform log(a^b) = b * log(a) to calculate log(3^3^3^3) = (3^3^3) log(3) = 7625597484987 * log(3), which works out to about 3638334640024.09968557 if you take logs base 10. You'll notice that the integer part of that gives you the number of digits, and if you take 10^0.09968557, you end up with 1.2580143 or so. Wolfram worked it out to a few more digits than I did, but this is pretty basic stuff with logarithms and not as expensive as computing 3^3^3^3 in integer arithmetic.

They also give the "last few digits" as 6100739387, but that's easily done using modular exponentiation: a recent version of Python will instantly return the same value for pow(3, 3**3**3, 10_000_000_000). Even though the power is rather large, the numbers being multiplied never get more than 10 digits long, so everything is easy to work out, and repeated squaring provides a major shortcut for the exponentiation.

Answer 2

This breaks both long long (integer, 18-19 digits) and double

the g++ compiler also provides a 128-bit type __int128_t with a value range of ?2^127 ... 2^127 ? 1 (about ?10^38 ... 10^38)

3^(7.625e+13) takes abnormally very long time

as you just calculate pow like this: return Int(base.val ** exp.val) guess, it takes O(N)

you could optimize it with fast powering algorithm ( or binary exponentiation )

    def bin_pow( a, n ):
        """
        calculates a ^ n
        complexity O(log n)
        a -> integer
        b -> integer
        ans -> integer
        """
        ans = 1
        while n:
          if( n % 2 != 0 ):
            ans *= a
          a*=a
          n /= 2
    return ans

or https://github.com/ampl/gsl is another alternative for C/C++

Answer 3

WolframAlpha shows the first few digits and the number of digits:

1.25801429062749131786039069820328121551804671431659... × 10^3638334640024

We can do that ourselves in less than a millisecond by keeping just the first let's say 50 digits. Do it once where we round down all the time, so we get a lower bound of the real value. Do it again where we round up all the time so we get an upper bound of the real value. The digits where lower and upper bound match are correct:

lower bound: 1.2580142906274913178603906982032812155096427774056 × 10^3638334640024
correct:     1.2580142906274913178603906982032812155... × 10^3638334640024
upper bound: 1.2580142906274913178603906982032812155218895229290 × 10^3638334640024
0.0003557720046956092 seconds

If you want more correct digits, just keep more digits and the lower and upper bound will match for more digits. For example with keep = 100:

lower bound: 1.258014290627491317860390698203281215518046714316596015189674944381211011300017785310803884345530895 × 10^3638334640024
correct:     1.258014290627491317860390698203281215518046714316596015189674944381211011300017785310803... × 10^3638334640024
upper bound: 1.258014290627491317860390698203281215518046714316596015189674944381211011300017785310803933426563879 × 10^3638334640024
0.0004895620222669095 seconds

My code (Try it online!) represents a number as two integers mantissa and exponent, representing the number mantissa × 10^exponent. Whenever mantissa grows too large (more than keep digits), it removes digits from mantissa and increases exponent accordingly:

def power(b, e):
    '''Compute b to the power of e.'''
    result = Int(1)
    while e:
        if e % 2:
            result *= b
        b *= b
        e //= 2
    return result

class Int:
    def __init__(self, mantissa, exponent=0):
        n = len(str(mantissa))
        if n > keep:
            remove = n - keep
            if round == 'down':
                mantissa //= 10**remove
            else:
                mantissa = -(-mantissa // 10**remove)
            exponent += remove
        self.mantissa = mantissa
        self.exponent = exponent
    def __mul__(self, other):
        return Int(self.mantissa * other.mantissa,
                   self.exponent + other.exponent)
    def __str__(self):
        m = str(self.mantissa)
        e = self.exponent
        if len(m) > 1:
            e += len(m) - 1
            m = m[0] + '.' + m[1:]
        return f'{m} × 10^{e}'

from os.path import commonprefix
from timeit import default_timer as timer
t0 = timer()

keep = 50
round = 'down'
lower = str(power(Int(3), 3**3**3))
round = 'up'
upper = str(power(Int(3), 3**3**3))

m, e = lower.split(' × ')
M, E = upper.split(' × ')
assert e == E
m = commonprefix([m, M])

print('lower bound:', lower)
print('correct:    ', m + '...', '×', e)
print('upper bound:', upper)

print(timer() - t0, 'seconds')