In the pass I did something similarly. At the time I applied Gauss-Legendre quadrature to evaluate $H_{n,r+1}$ intergals; this calculation adopted (1) Gauss-Legendre weight and points calculated in advance (then, change the variable $t$ in the integrand to $\xi \in [-1,1]$), and (2) multi-precision floating point library fmlib.
I guess that you obtain $H_{n,i} $ from Lagrange Interpolation Formulae. Let me review some basic thing so that you may recognize the same derivation. From the Lagrange Interpolation, we have
\begin{align}
f(x)
=
\sum_{i=1}^{n}
l_i (x)
f_i
+
R_n (x)
\end{align}
where
\begin{align}
l_i (x)
&=
\frac{
\Pi_n( x)
}{
(x-x_i)
}
.
\frac{
1
}{
\Pi _n ' (x_i)
}
=
\frac{
( x - x_1 )
\ldots
( x - x_{i-1} )
( x - x_{i+1} )
\ldots
( x - x_{n} )
}{
( x_i - x_1 )
\ldots
( x_i - x_{i-1} )
( x_i - x_{i+1} )
\ldots
( x_i - x_{n} )
} ,
\\
R_n (x)
&=
\frac{ f^{(n)} (\xi) }{ n! } \cdot
\Pi _n (x)
,
\end{align}
for some $\xi = \xi(x) \in (x_1, x_n )$.
Note that
we defined
$
\Pi_n (x)
=
\prod_{i=1}^{n}
(x-x_i)
$,
and we also have
\begin{align}
\left| R_n (x) \right|
&\le
\frac{ ( x_n - x_1 )^{n} }{ n! }
\max_{ x_1 \le x \le x_n } \left| f^{(n)} (x) \right|
\end{align}
which can be proved using the Role’s theorem.
Now we take
\begin{align}
\int_{x_1}^{x_n}
f(x)
d x
&=
\sum_{i=1}^{n}
\underbrace{
\left(
\int_{x_1}^{x_n}
l_i (x)
d x
\right)
}_{ A_ i }
f_i
+
\underbrace{
\int_{x_1}^{x_n}
R_n (x)
d x
}_{ B }
\end{align}
Define
\begin{align}
A_ i
&=
\int_{x_1}^{x_n}
l_i (x)
d x ~~~(\equiv H_{n,i} ? )
\\
B
&=
\int_{x_1}^{x_n}
R_n (x)
d x
\end{align}
We are going to estimate those integrals. Change variable,
\begin{align}
h
&=
\frac{ x_n - x_1 }{ n - 1 }
,~~~~~
x = x_1 + (t-1) h
,~~
\forall t \in [1,n]
\end{align}
Using the latter one, $x - x_j = x_1 + (t-1)h - x_j $ and
$ x_j = x_1 + (j-1) h $,
we have
$$
x - x_j
=
x_1 + (t-1)h -
\left[ x_1 + (j-1)h \right]
=
(t-j)h
$$
Thus
\begin{align}
A_ i
&=
\int_{x_1}^{x_n}
l_i (x)
d x
=
\int_{x_1}^{x_n}
\frac{
( x - x_1 )
\ldots
( x - x_{i-1} )
( x - x_{i+1} )
\ldots
( x - x_{n} )
}{
( x_i - x_1 )
\ldots
( x_i - x_{i-1} )
( x_i - x_{i+1} )
\ldots
( x_i - x_{n} )
}
d x
\\
&=
\int_{1}^{n}
\frac{
h ( t - 1 )
\cdot
h ( t - 2 )
\cdots
h ( t-i+1 )
\cdot
h (t-i-1)
\cdots
h(t-n)
}{
h ( i - 1 )
\cdot
h ( i - 2 )
\cdots
h ( i-i+1 )
\cdot
h (i-i-1)
\cdots
h(i-n)
}
\cdot
h d t
\end{align}
Or
\begin{align}
A_ i
&=
h
\int_{1}^{n}
\frac{
( t - 1 )
\cdot
( t - 2 )
\cdots
( t-i+1 )
\cdot
(t-i-1)
\cdots
(t-n)
}{
( i - 1 )
\cdot
( i - 2 )
\cdots
( i-i+1 )
\cdot
(i-i-1)
\cdots
(i-n)
}
d t
\end{align}
And
\begin{align}
B
&=
\int_{x_1}^{x_n}
R_n (x)
d x
=
\int_{x_1}^{x_n}
( x - x_1 )
\ldots
( x - x_{i} )
\ldots
( x - x_{n} )
\cdot
\frac{ f^{(n)} (\xi (x) ) }{ n! }
d x
\\
&=
\int_{1}^{n}
h ( t - 1 )
\cdot
h ( t - 2 )
\cdots
h ( t - i)
\cdots
h(t-n)
\cdot
\frac{ f^{(n)} (\xi (x(t)) ) }{ n! }
\cdot
h d t
\end{align}
We estimate
\begin{align}
| B |
&\le
M
h^{n+1}
\end{align}
where $M = |C_n | \sup_{ \xi \in [x_1,x_n] } \left| f^{(n)} (\xi) \right| $
and
\begin{align}
C_n
&=
\frac{1}{ n! }
\int_{1}^{n}
( t - 1 )
\cdots
( t - i)
\cdots
( t - n )
d t
\label{c}
\end{align}
Here integrals $A_i$ and $C_n$ can be estimated EXACTLY using
Gauss-Legendre formulas with $n/2$ points.
As I could see in numerical result at the time, we have (a) $A_i $ are increased to very large values as $n$ increasing, for example $\max A_i \sim 10^{0}, 10^{2}, 10^{4}, 10^{7}, 10^{10}$ and $10^{24}$ when $n=10,20,30,40,50$, and $100$, respectively; (b) $A_i$ is symmetry with repect to $i$, i.e. $A_{i} = A_{n-i+1}$; (c) $\sum_{i=1}^{n} A_{i} = n-1$; and (c) $C_n = 0$ for $n$ odd, and $C_n$ is bounded by a small value when $n$ even.
