How to reasonably (numerically) estimate $n\int_0^\infty\left(1-\left(1-e^{-t}\left(1+t+\ldots+{{t^{m-1}}\over{(m-1)!}}\right)\right)^n\right)dt$? Recently in doing some expected value calculations, I've derived the following two integrals:
$$6\int_0^\infty\left(1-\left(1-e^{-t}\left(1+t\right)\right)^6\right)dt$$
$$6\int_0^\infty\left(1-\left(1-e^{-t}\left(1+t + {{t^2}\over2}\right)\right)^6\right)dt$$
Plugging these into Wolfram Alpha can get us numerical answers. For instance, the first integral comes out to about $24.1$, and the second integral comes out to about $32.7$.
However, I'm wondering if anyone can give a reasonable numerical estimate for these $2$ integrals from first principles (pencil and paper) without using a calculator, Wolfram Alpha, or a computer.
I've tried but made little to no progress, and I consulted some nearby PhD students and they didn't know either, so asking here.

I doubt that you can get a reasonable estimate with pencil and paper.
For a quick and rough approximation (still requires calculator):
Let $E_1 = 6 I_1 = 6 \int_{0}^\infty g_1(t) dt$ with $g_1(t)=1-h_1(t)^6$ and $h_1(t)=1-e^{-t}(1+t)$.
Same for $E_2$, with $h_2(t)=1-e^{-t}(1+t+\frac12 t^2)$
In both cases, we can see that $h(t)$ increase monotonously from $0$ to $1$ (can be checked computing the derivative). Hence we can expect that, after the composition with $1-(\cdot)^6$, $g(t)$ has a sigmoid shape, starting at $g(0) =1$ with some quite steep transition to $g(+\infty)=0$. If we assume this, we can try to approximate this sigmoid by a piecewise linear activation function. And in this case, the integral is $I=t_{0}$, where $g(t_0)=1/2$.
Under this assumption, we are left with computing $t_0$ as the solution of
$$  1-(1-e^{-t}(1+t))^6=1/2 \tag 1$$
or $$ t = \log(1+t) - \log \left(1- \sqrt[6]{1/2}\right) \tag 2$$
Again, this cannot be computed with paper and pencil, but the numerical solution can be found iteratively, eg with a spreadsheet.
(The attempt to use the Taylor series for the logarithm leads nowhere).
We get the values $I_1=3.78$ ($E_1=22.68$) for the first one
and $I_2=5.196$ ($E_2=31.17$) for the second one.
The relative error wrt the exact values (from Yuri Negometyanov's answer) is around $5\%$.
See also  Qiaochu Yuan's comment for a probabilistic bound.
For general $m,n$, the method used in the answer of this similar (but not quite) question might be of interest.
