## Calculating Sophomore’s dream

**Posted: 2020-07-16 12:00**, last updated: 2022-03-02 20:39.

Sophomore's dream is a mathematical constant defined by

$$ S = \sum_{n = 1}^{\infty} n^{-n}. $$

This constant converges to 1.291285997…

I hold the (unaccredited) world record for computing the most decimal digits of $S$. As of September 2017, I have computed 1,000,000 decimal digits of $S$, thereby dwarfing my previous world record of 200,000 decimal digits (2013).

Here is a frequency table:

digit | 0 | 1 | 2 | 3 | 4 |
---|---|---|---|---|---|

# of occurences | 100194 | 100295 | 100037 | 99783 | 100162 |

digit | 5 | 6 | 7 | 8 | 9 |

# of occurences | 99949 | 100054 | 99937 | 99855 | 99734 |

The constant is called *Sophomore's dream* because early students of mathematics, when first exposed to the constant, are dismayed at the fact that there is neither an easy algorithm to compute it nor any closed form, despite its pleasant and straightforward definition. In fact, it is one of the most difficult to compute commonly known mathematical constants.

I want to share how I computed so many digits and what are some pitfalls.

First of all, I compute the individual terms. This is easier said than done. Importantly, it is not possible to use builtin integer types (at least not those of any of my preferred programming languages) because they lack precision. Think of $712583^{712583}$: This number is too large to be expressed precisely by `int`

. Even the largest unsigned integers currently implemented in orthodox C++, `uint64_t`

, cannot hold these large values. `float`

and others fail equally (and they should not be used here regardless). Therefore, I used GMP to achieve arbitrary precision.

I then calculated the inverse of all the terms, i. e. $n^n$ instead of $n^{-n}$. I did this using a custom `pow()`

function.

However, it is of course not possible to calculate *all* the terms since there are infinitely many. I had to figure out how many terms I had to calculate to get at least 1,000,000 decimal digits. In other words, I had to find the solution to the equation

$$ \log_{10}\left(x^x\right) = 1\,000\,000. $$

The solution is $x = 189481.28$ terms. To be safe, I "rounded" that up to $x = 190500$, which should actually give me 1,005,820 digits – again, to be safe. This calculation is allowed since the constant converges rapidly.

Thereafter, I created a custom GMP implementation of long division to divide 1 by the inverse terms, trivially giving me the actual terms $n^{-n}$ with a specified precision.

Then, I added all the terms that I computed.

By the way: Depending on the number of terms, it may be prudent to use my p-add application to add the individual fractions (i.e. $\frac{1}{n^n}$). I used this method when I calculated 200,000 decimal digits. You repeat until you are left with just one very large fraction that is then expressed numerically using long division. You essentially exponentiate, add the fractions and divide. However, using this method is not a good idea for a much larger number of decimal digits since it quickly requires prohibitive amounts of continuous random access memory. You then have to refer to the above methodology (i.e. exponentiate, divide, add).

Sounds easy? Well, it sort of is. However, this easy process involves developing many highly performant, multi-threaded C++ tools. Most important, however, is the time and space constraint: I estimate that my computation took me, all in all, about 240 hours and more than 300 GiB of disk space. On the other hand, using the *exponentiate, divide, add* methodology does not require immense amounts of memory. It is a nice project for the time off and I therefore highly encourage you to surpass my world record.