Computing a definite integral for the exponential of a polynomial. I'm looking for the most effective way to compute the following definite integral for a given n:
$\int_0^1e^{2\sum_{i=1}^n\frac{x^i-1}{i}}dx$
There are elementary solutions for n = 0 and n = 1, but I'm aware larger values of n can only be solved with non-elementary functions. For example, there is a solution for n = 2 using the imaginary error function.
My questions are:

*

*Are there any non-elementary functions that could provide a generic solution for an arbitrary n.

*Alternatively, is it perhaps possible to solve the integral in the form of a Taylor Series?

*Alternatively, are there any tricks that could help me compute a definite integral faster for instances of this case?

In short, I'm trying to see what options I have for computing the above integral for arbitrary n, short of just using standard tools for computing the value of a definite integral.
 A: Numerical integration
According to the Gauss–Kronrod quadrature method:
$$
I\left(k_{max}\right) := \int_{-1}^1 \underbrace{\exp\left[2\sum_{k=1}^{k_{max}}\frac{\left(\frac{\tau + 1}{2}\right)^k - 1}{k}\right]\frac{1}{2}}_{:= f(\tau)}\,\text{d}\tau \approx \sum_{i = 1}^n f(t_i)\,w_i
$$
where the respective nodes and weights can be tabulated directly in Wolfram Mathematica 12.2 as indicated in detail on this page, or in other spreadsheets such as Matlab in this other page.
As an example, fixed $k_{max} = 20$, already with the classic:
f[τ_] := Exp[2 Total[Table[(((τ + 1)/2)^k - 1)/k, {k, 20}]]] / 2
{n, prec} = {7, 10};

(*Gauss-Legendre quadrature*)
t = 2 NIntegrate`GaussRuleData[n, prec][[1]] - 1;
w = 2 NIntegrate`GaussRuleData[n, prec][[2]];
IG = Total[Table[f[t[[i]]] w[[i]], {i, n}]]

(*Gauss-Kronrod quadrature*)
t = 2 NIntegrate`GaussKronrodRuleData[n, prec][[1]] - 1;
w = 2 NIntegrate`GaussKronrodRuleData[n, prec][[2]];
IK = Total[Table[f[t[[i]]] w[[i]], {i, 2 n + 1}]]

Abs[IG - IK]

$$
I_{G7} = 0.034964534
\; \; \; \; \; \;
I_{K15} = 0.036094679
\; \; \; \; \; \;
\left|I_{G7} - I_{K15}\right| = 1.1 \cdot 10^{-3}
$$
the percentage error is around $1.6\%$, not bad! On the other hand, slightly increasing the nodes:
$$
I_{G15} = 0.036093780
\; \; \; \; \; \;
I_{K31} = 0.036094015
\; \; \; \; \; \;
\left|I_{G15} - I_{K31}\right| = 2.4 \cdot 10^{-7}
$$
the percentage error is practically zero and therefore $I \approx I_{K31}$, which is what is required.

Curve fitting
Keeping the number of nodes and the calculation precision unchanged, writing:
t = 2 NIntegrate`GaussKronrodRuleData[15, 10][[1]] - 1;
w = 2 NIntegrate`GaussKronrodRuleData[15, 10][[2]];

data1 = Table[f[τ_] := Exp[2 Total[Table[(((τ + 1)/2)^k - 1)/k, {k, kmax}]]] / 2;
              Total[Table[f[t[[i]]] w[[i]], {i, 31}]], {kmax, 0, 100}];

data2 = Table[0.74/(kmax + 1) + 0.26/(kmax + 1)^2, {kmax, 0, 100}];

ListPlot[{data1, data2}, AxesLabel -> {"kmax", "I(kmax)"}, PlotStyle -> {Blue, Red}]

we get:

from which it's clear that function estimates the integral very well (maximum error $0.5\%$).
