Roll to see all sides!
Python + NumPy, 3.95 Seconds



from __future__ import division
import numpy as np

def rolls(n):
    if n == 1:
        return 1
    r = n * (np.log(n) - np.log(np.log(2)))
    x = np.log1p(np.arange(n) / -n)
    cx = x.cumsum()
    y = cx[:-1] + cx[-2::-1] - cx[-1]
    while True:
        r0 = np.round(r)
        z = np.exp(y + r0 * x[1:])
        z[::2] *= -1
        r = r0 - (z.sum() + 0.5) / z.dot(x[1:])
        if abs(r - r0) &lt; 0.75:
            return np.ceil(r).astype(int)

for n in [1, 2, 3, 4, 5, 6, 20, 52, 100, 366, 1000, 2000, 5000, 10000, 100000, 1000000]:
    print('R({}) = {}'.format(n, rolls(n)))

import timeit
print('Benchmark: {:.2f}s'.format(timeit.timeit(lambda: sum(map(rolls, range(1, 10001))), number=1)))


Try it online!

How it works

This uses the closed-form series for P(r, n), and its derivative with respect to r, rearranged for numerical stability and vectorization, to do a Newton’s method search for r such that P(r, n) = 0.5, rounding r to an integer before each step, until the step moves r by less than 3/4.  With a good initial guess, this usually takes just one or two iterations.

xi = log(1 − i/n) = log((n − i)/n)
cxi = log(n!/((n − i − 1)! ⋅ ni + 1)
yi = cxi + cxn − i − 2 − cxn − 1 = log binom(n, i + 1)
zi = (-1)i + 1 ⋅ binom(n, i + 1) ⋅ ((n − i − 1)/n)r
1 + ∑zi = n! ⋅ S(r, n)/nr = P(r, n)
zi ⋅ xi + 1 = (-1)i + 1 ⋅ binom(n, i + 1) ⋅ ((n − i − 1)/n)r log((n − i − 1)/n)
∑zi ⋅ xi + 1 = d/dr P(r, n)