There is a simple closed formula in terms of the roots of a degree-6 polynomial.
It's actually a little easier to consider a general fair die with $d\ge 2$ faces labeled with the numbers $1,2,\ldots, d.$
Let $e_k$ be the expected number of rolls needed to equal or exceed $k.$ For $k\le 0,$ $e_k=0.$ Otherwise the expectation is one more than the expectation of the number of rolls to reach the immediately preceding value, which would be among $k-d, k-d+1, \ldots, k-1,$ whence
$$e_k = 1 + \frac{1}{d}\left(e_{k-d} + e_{k-d+1} + \cdots + e_{k-1}\right).\tag{1}$$
This linear recurrence relation has a solution in the form
$$e_k = \frac{2k}{d+1} + \sum_{i=1}^d a_i \lambda_i^k\tag{2}$$
where the $\lambda_i$ are the $d$ complex roots of the polynomial
$$T^d - \frac{1}{d}(T^{d-1} + T^{d-2} + \cdots + T + 1).\tag{3}$$
The constants $a_i$ are found by applying the solution $(2)$ to the values $k=-(d-1), -(d-2), \ldots, -1, 0$ where $e_k=0$ in every case. This gives a set of $d$ linear equations in the $d$ constants and it has a unique solution. That the solution works can be demonstrated by verifying the recurrence $(1)$ using the fact that every root satisfies $(3):$
$$\eqalign{
1 + \frac{1}{d}\sum_{j=1}^{d} e_{k-j} &= 1 + \frac{1}{d}\sum_{j=1}^{d} \left(\frac{2(k-j)}{d+1} + \sum_{i=1}^d a_i \lambda_i^{k-j}\right) \\
&= \frac{2k}{d+1} + \sum_{i=1}^d a_i \lambda_i^{k-d}\left[\frac{1}{d}(1 + \lambda_i + \cdots + \lambda_i^{d-1})\right] \\
&= \frac{2k}{d+1} + \sum_{i=1}^d a_i \lambda_i^{k-d}\lambda_i^d \\
&= \frac{2k}{d+1} + \sum_{i=1}^d a_i \lambda_i^k = e_k.
}$$
This closed form solution gives us good ways to approximate the answer as well as to evaluate it accurately. (For small to modest values of $k,$ direct application of the recurrence is an effective computational technique.) For example, with $d=6$ we can readily compute
$$e_{1\,000\,000} = 285714.761905\ldots$$
For approximations, there will be a unique largest root $\lambda_{+}=1$ so eventually (for sufficiently large $k$) the term $\lambda_{+}^k$ will dominate the $d$ terms in $(2).$ The error will decrease exponentially according to the second smallest norm of the roots. Continuing the example with $k=6,$ the coefficient of $\lambda_{+}$ is $a_{+}=0.4761905$ and the next-smallest norm is $0.7302500.$ (Incidentally, the other $a_i$ tend to be very close to $1$ in size.) Thus we may approximate the previous value as
$$e_{1\,000\,000} \approx \frac{2\times 10^6}{6+1} + 0.4761905 = 285714.761905\ldots$$
with an error on the order of $0.7302500^{10^6} \approx 10^{-314\,368}.$
To demonstrate how practical this solution is, here is R code that returns a function to evaluate $e_k$ for any $k$ (within the scope of double precision floating point calculations) and not overly large $d$ (it will bog down once $d\gg 100$):
die <- function(d, mult=1, cnst=1, start=rep(0,d)) {
# Create the companion matrix (its eigenvalues are the lambdas).
X <- matrix(c(0,1,rep(0,d-1)),d,d+1)
X[, d] <- mult/d
lambda <- eigen(X[, 1:d], symmetric=FALSE, only.values=TRUE)$values
# Find the coefficients that agree with the starting values.
u <- 2*cnst/(d+1)
a <- solve(t(outer(lambda, 1:d, `^`)), start - u*((1-d):0))
# This function assumes the starting values are all real numbers.
f <- Vectorize(function(i) Re(sum(a * lambda ^ (i+d))) + u*i)
list(f=f, lambda=lambda, a=a, multiplier=mult, offset=cnst)
}
As an example of its use, here it calculates the expectations for $k=1,2,\ldots, 16:$
round(die(6)$f(1:10), 3)
1.000 1.167 1.361 1.588 1.853 2.161 2.522 2.775 3.043 3.324 3.613 3.906 4.197 4.476 4.760 5.046
The object it returns includes the roots $\lambda_i$ and their multipliers $a_i$ for further analysis. The first component of the multipliers array is the useful coefficient $a_{+}.$
(If you're curious what the other parameters of die are for, execute die(2, 2, 0, c(1,0))$f(1:10) and see whether you recognize the output ;-). This generalization assisted in developing and testing the function.)
