# The number of possible numeric outcomes of parenthesizations of 2^2^…^2

Consider an expression 2^2^...^2 with n operators ^. Operator ^ means exponentiation ("to the power of"). Assume that it has no default assosiativity, so the expression needs to be fully parenthesized to become unambiguous. The number of ways to parenthesize the expression are given by Catalan numbers C_n=(2n)!/(n+1)!/n!.

Sometimes different parenthesizations give the same numeric outcome, for example (2^2)^(2^2)=((2^2)^2)^2, so the number of different possible numeric outcomes for a given n is less than C_n for all n>1. The sequence starts as 1, 1, 2, 4, 8, ... as opposed to Catalan numbers 1, 2, 5, 14, 42, ...

The problem is to write the fastest program (or function) that accepts n as an input and returns the number of different possible numeric outcomes of the expression 2^2^...^2 with n operators ^. The performance should not significantly deteriorate as n grows, so direct calculation of high power towers is probably a bad idea.

• I'm just sharing an idea here, but it seems that it should be possible to exclusively use addition and multiplication, as the answer will always be of the form 2^n, and it therefore would be unnecessary to keep track of anything except n. I.e., just using the rules of exponentiation seems wise. However, there surely is a smarter and completely algebraic way to do this. – Fors May 9 '13 at 15:36
• @Fors I guess n is still way too big to compute. Still, well noted. Maybe a recursive representation in the form "1 or 2^(...) or (...)+(...)"; but you still have the problem of how to normalize such representation of a number (or compare two representations for value equality). – John Dvorak May 9 '13 at 16:24
• @JanDvorak, A002845 (no closed form given) – Peter Taylor May 9 '13 at 17:44
• Related oeis.org/A003018/a003018.pdf – Dr. belisarius May 9 '13 at 18:00
• @Vladimir Reshetnikov: I think there is an off-by-one error in your formula. When you have n twos and C_n=(2n)!/(n+1)!/n! should be the number of parenthesizations, then for n=3 it should be 5, correct? I see (2^2)^2 and 2^(2^2), but what are the other three combinations? I think C_n gives you the number of parenthesizations for n+1 twos. – Martin Thoma May 15 '13 at 5:37

# Python 2.7

This approach takes advantage of the following considerations:

Any integer can be represented as a sum of powers of two. The exponents in the powers of two can also be represented as powers of two. For example:

8 = 2^3 = 2^(2^1 + 2^0) = 2^(2^(2^0) + 2^0)


These expressions that we end up with can be represented as sets of sets (in Python, I used the built-in frozenset):

• 0 becomes the empty set {}.
• 2^a becomes the set containing the set representing a. E.g.: 1 = 2^0 -> {{}} and 2 = 2^(2^0) -> {{{}}}.
• a+b becomes the concatenation of the sets representing a and b. E.g., 3 = 2^(2^0) + 2^0 -> {{{}},{}}

It turns out that expressions of the form 2^2^...^2 can easily be transformed into their unique set representation, even when the numerical value is much too large to be stored as an integer.

For n=20, this runs in 8.7s on CPython 2.7.5 on my machine (a bit slower in Python 3 and much slower in PyPy):

"""Analyze the expressions given by parenthesizations of 2^2^...^2.

Set representation:  s is a set of sets which represents an integer n.  n is
given by the sum of all 2^m for the numbers m represented by the sets
contained in s.  The empty set stands for the value 0.  Each number has
exactly one set representation.

In Python, frozensets are used for set representation.

Definition in Python code:
def numeric_value(s):
n = sum(2**numeric_value(t) for t in s)
return n"""

import itertools

def single_arg_memoize(func):
"""Fast memoization decorator for a function taking a single argument.

The metadata of <func> is *not* preserved."""

class Cache(dict):
def __missing__(self, key):
self[key] = result = func(key)
return result
return Cache().__getitem__

def count_results(num_exponentiations):
"""Return the number of results given by parenthesizations of 2^2^...^2."""
return len(get_results(num_exponentiations))

@single_arg_memoize
def get_results(num_exponentiations):
"""Return a set of all results given by parenthesizations of 2^2^...^2.

<num_exponentiations> is the number of exponentiation operators in the
parenthesized expressions.

The result of each parenthesized expression is given as a set.  The
expression evaluates to 2^(2^n), where n is the number represented by the
given set in set representation."""

# The result of the expression "2" (0 exponentiations) is represented by
# the empty set, since 2 = 2^(2^0).
if num_exponentiations == 0:
return {frozenset()}

# Split the expression 2^2^...^2 at each of the first half of
# exponentiation operators and parenthesize each side of the expession.
split_points = xrange(num_exponentiations)
splits = itertools.izip(split_points, reversed(split_points))
splits_half = ((left_part, right_part) for left_part, right_part in splits
if left_part <= right_part)

results = set()
for left_part, right_part in splits_half:
for left in get_results(left_part):
for right in get_results(right_part):
return results

def exponentiate(base, exponent):
"""Return the result of the exponentiation of <operands>.

<operands> is a tuple of <base> and <exponent>.  The operators are each
given as the set representation of n, where 2^(2^n) is the value the
operator stands for.

The return value is the set representation of r, where 2^(2^r) is the
result of the exponentiation."""

# Where b is the number represented by <base>, e is the number represented
# by <exponent> and r is the number represented by the return value:
#   2^(2^r) = (2^(2^b)) ^ (2^(2^e))
#   2^(2^r) = 2^(2^b * 2^(2^e))
#   2^(2^r) = 2^(2^(b + 2^e))
#   r = b + 2^e

# If <exponent> is not in <base>, insert it to arrive at the set with the
# value: b + 2^e.  If <exponent> is already in <base>, take it out,
# increment e by 1 and repeat from the start to eventually arrive at:
#   b - 2^e + 2^(e+1) =
#   b + 2^e
while exponent in base:
base -= {exponent}
exponent = successor(exponent)
return base | {exponent}

@single_arg_memoize
def successor(value):
"""Return the successor of <value> in set representation."""
# Call exponentiate() with <value> as base and the empty set as exponent to
# get the set representing (n being the number represented by <value>):
#   n + 2^0
#   n + 1
return exponentiate(value, frozenset())

def main():
import timeit
print timeit.timeit(lambda: count_results(20), number=1)
for i in xrange(21):
print '{:.<2}..{:.>9}'.format(i, count_results(i))

if __name__ == '__main__':
main()


(The memoization decorator's concept is copied from http://code.activestate.com/recipes/578231-probably-the-fastest-memoization-decorator-in-the-/ .)

Output:

8.667753234
0...........1
1...........1
2...........1
3...........2
4...........4
5...........8
6..........17
[...]
19.....688366
20....1619087


Timings for different n:

 n    time
16    0.240
17    0.592
18    1.426
19    3.559
20    8.668
21   21.402


Any n above 21 results in a memory error on my machine.

I'd be interested if anyone can make this faster by translating it into a different language.

Edit: Optimized the get_results function. Also, using Python 2.7.5 instead of 2.7.2 made it run a bit faster.

• I made a C# translation but using sorted arrays and doing the addition in order rather than by set contains checks. It's much slower, and I haven't yet profiled to see whether that's due to not memoising the successor function or due to the cost of comparisons. – Peter Taylor Jul 24 '13 at 17:42
• I haven't profiled @flornquake's (brilliant) code, but I'd assume much of the CPU time is spent doing set membership tests and set manipulation operations, which are both fairly well optimized in Python, using its ubiquitous hash table and hash key routines. Memoization is certainly a big thing, with an exponential algorithm such as this. If you left that out, you can expect exponentially slower performance. – Tobia Jul 27 '13 at 22:39
• @Tobia, actually I found that in C# memoising the successor function made it slower. I also found that a more literal translation (using set operations) was significantly slower than my lower level addition. The only real improvement I've found over my original code was to take into account (a^b)^c = (a^c)^b, and it's still much slower than this Python implementation. – Peter Taylor Jul 29 '13 at 11:02
• @PeterTaylor: Edit: As far as I can see, flornquake's algorithm relies on building sets of trees, where a tree is a set of trees itself, and so on. All the pieces of these trees, from the smallest empty set to the largest set of sets, are memoized. This means that all these trees contain "repeated structure" that is only computed once (by the CPU) and stored once (in RAM). Are you sure that your "addition in order" algorithm is identifying all of this repeated structure and computing it once? (what I called exponential complexity above) See also en.wikipedia.org/wiki/Dynamic_programming – Tobia Jul 29 '13 at 15:16
• @Tobia, we overlapped. I've posted the code. – Peter Taylor Jul 29 '13 at 15:18

## C#

This is a translation of flornquake's Python code to C# using a lower level addition routine which provides a moderate speedup over a direct translation. It's not the most optimised version I have, but that is quite a bit longer because it has to store the tree structure as well as the values.

using System;
using System.Collections.Generic;
using System.Linq;

namespace Sandbox {
class PowerTowers {
public static void Main() {
DateTime start = DateTime.UtcNow;
for (int i = 0; i < 17; i++)
Console.WriteLine("{2}: {0} (in {1})", Results(i).Count, DateTime.UtcNow - start, i);
}

private static IList<HashSet<Number>> _MemoisedResults;

static HashSet<Number> Results(int numExponentations) {
if (_MemoisedResults == null) {
_MemoisedResults = new List<HashSet<Number>>();
_MemoisedResults.Add(new HashSet<Number>(new Number[] { Number.Zero }));
}

if (numExponentations < _MemoisedResults.Count) return _MemoisedResults[numExponentations];

HashSet<Number> rv = new HashSet<Number>();
for (int i = 0; i < numExponentations; i++) {
IEnumerable<Number> rhs = Results(numExponentations - 1 - i);
foreach (var b in Results(i))
foreach (var e in rhs) {
}
}
return rv;
}
}

// Immutable
struct Number : IComparable<Number> {
public static Number Zero = new Number(new Number);
public static Number One = new Number(Zero);

// Ascending order

private Number(params Number[] children) {
_Children = children;
_Depth = children.Length == 0 ? 0 : 1 + children[children.Length - 1]._Depth;

int hashCode = 0;
foreach (var n in _Children) hashCode = hashCode * 37 + n.GetHashCode() + 1;
_HashCode = hashCode;
}

// Work forwards because children are in ascending order.
int off1 = 0, off2 = 0;
IList<Number> result = new List<Number>();
Number? carry = default(Number?);

while (true) {
if (!carry.HasValue) {
// Simple case
if (off1 < _Children.Length) {
if (off2 < n._Children.Length) {
int cmp = _Children[off1].CompareTo(n._Children[off2]);
else if (cmp == 0) {
off2++;
}
}
}
else if (off2 < n._Children.Length) result.Add(n._Children[off2++]);
else return new Number(result.ToArray()); // nothing left to add
}
else {
// carry is the (possibly joint) smallest value
int matches = 0;
if (off1 < _Children.Length && carry.Value.Equals(_Children[off1])) {
matches++;
off1++;
}
if (off2 < n._Children.Length && carry.Value.Equals(n._Children[off2])) {
matches++;
off2++;
}

if ((matches & 1) == 0) result.Add(carry.Value);
carry = matches == 0 ? default(Number?) : carry.Value.Add(One);
}
}
}

public Number Exp2() {
return new Number(this);
}

public int CompareTo(Number other) {
if (_Depth != other._Depth) return _Depth.CompareTo(other._Depth);

// Work backwards because children are in ascending order
int off1 = _Children.Length - 1, off2 = other._Children.Length - 1;
while (off1 >= 0 && off2 >= 0) {
int cmp = _Children[off1--].CompareTo(other._Children[off2--]);
if (cmp != 0) return cmp;
}

return off1.CompareTo(off2);
}

public override bool Equals(object obj) {
if (!(obj is Number)) return false;

Number n = (Number)obj;
if (n._HashCode != _HashCode || n._Depth != _Depth || n._Children.Length != _Children.Length) return false;
for (int i = 0; i < _Children.Length; i++) {
if (!_Children[i].Equals(n._Children[i])) return false;
}

return true;
}

public override int GetHashCode() {
return _HashCode;
}
}
}