A question about generalized harmonic numbers modulo $p$ Let $p \equiv 1 \pmod{3}$ be a prime and denote $H_{n,m} = \sum_{k = 1}^n 1/k^m$ as the $n,m$-th generalized harmonic number.  I'm interested in computing $H_{(p-1)/3,\, 2}$ and $H_{(p-1)/6,\,2}$ modulo $p$.  From this paper I know
\begin{align*}
H_{(p-1)/6,2} \equiv -\frac{B_{2p-3}(1/6)}{2p-3}\pmod p
\end{align*}
and
\begin{align*}
H_{(p-1)/3,2} \equiv -\frac{B_{2p-3}(1/3)}{2p-3}\pmod p
\end{align*}
where $B_n(x) = \sum_{k = 0}^n {n \choose k}x^{n-k}B_k$ is the $n$-th Bernoulli polynomial.
It is well known, I think, that
\begin{align*}
B_n(1/6) = (1 - 3^{1-n})(1 - 2^{1-n})\frac{B_n}{2\cdot 3^{n-1}}
\end{align*}
and
\begin{align*}
B_n(1/3)=(1 - 3^{1-n})\frac{B_n}{2\cdot 3^{n-1}}
\end{align*}
but only for even $n$.  From this we can eventually express $H_{(p-1)/3} = H_{(p-1)/3, 1}$ and $H_{(p-1)/6} = H_{(p-1)/6,1}$ modulo $p$ in terms of the Fermat quotients $q_p(2)$ and $q_p(3)$ where $q_p(a) = (a^{p-1} - 1)/p$ for $a$ co-prime to $p$.
It is mentioned in the paper above that there is some kind of expression for for $B_n(1/3)$ and $B_n(1/6)$ when $n$ is odd in terms of "$I$ numbers", but I can't find the source, or anything else that is relevant.
Is it possible to write $H_{(p-1)/3,2}$ and $H_{(p-1)/6,2}$ in terms of Fermat quotients modulo $p$ analogous to $H_{(p-1)/3}$ and $H_{(p-1)/6}$?  What are these "$I$ numbers?"?
I've asked a similar question here but I thought I would ask here.
 A: Glaisher's I-numbers are described in J. W. L. Glaisher, On a set of coefficients analogous to the Eulerian numbers, Proc. London Math. Soc., 31 (1899), 216-235.
A: Here is what I've learned:
\begin{align*}
2\sum_{n = 0}^{\infty} B_{2n+1}\left(\frac{1}{3}\right) \frac{x^{2n+1}}{(2n+1)!} & = \frac{xe^{\frac{1}{3}x}}{e^x -1} + \frac{xe^{-\frac{1}{3}x}}{e^{-x} -1} =  \frac{xe^{\frac{1}{3}x}}{e^x -1} - \frac{xe^{\frac{2}{3}x}}{e^{x} -1} 
\\
&= \frac{x(e^{\frac{1}{3}x}-e^{\frac{2}{3}x})}{e^x -1} = \frac{xe^{\frac{1}{3}x}(1-e^{\frac{1}{3}x})}{(e^{\frac{1}{3}x} -1)(e^{\frac{2}{3}x}+e^{\frac{1}{3}x}+1)}
\\
& = \frac{-xe^{\frac{1}{3}x}}{e^{\frac{2}{3}x}+e^{\frac{1}{3}x}+1} = \frac{-x}{e^{\frac{1}{3}x}+e^{-\frac{1}{3}x}+1} 
\end{align*}
The $n$-th $I$ number is defined by
\begin{align*}
\frac{3/2}{e^x + e^{-x} + 1} = \sum_{n = 0}^{\infty} I_n \frac{x^n}{n!}
\end{align*}
Notice that
\begin{align*}
2 \sum_{n = 0}^{\infty} I_{2n+1} \frac{x^{2n+1}}{(2n+1)!} &=  \sum_{n = 0}^{\infty} I_n \frac{x^n}{n!} -  \sum_{n = 0}^{\infty} I_n \frac{(-x)^n}{n!}
 \\
 & = 0.
 \end{align*}
Therefore, $I_n = 0$ whenever $n$ is odd, and $ \sum_{n = 0}^{\infty} I_n \frac{x^n}{n!} =  \sum_{n = 0}^{\infty} I_{2n} \frac{x^{2n}}{(2n)!}$.  Now,
\begin{align*}
 \frac{-x}{e^{\frac{1}{3}x}+e^{-\frac{1}{3}x}+1} &= -\frac{2}{3}x \frac{3/2}{e^{\frac{1}{3}x}+e^{-\frac{1}{3}x}+1}
   \\
  & = -\frac{2}{3}x \sum_{n = 0}^{\infty} I_{2n} \frac{\left( \frac{1}{3}x\right)^{2n}}{(2n)!}
  \\
  & = 2\sum_{n = 0}^{\infty} -\frac{(2n+1)I_{2n}}{3^{2n+1}} \frac{x^{2n+1}}{(2n+1)!}.
 \end{align*}
Therefore,
\begin{align*}
 \sum_{n = 0}^{\infty} -\frac{(2n+1)I_{2n}}{3^{2n+1}} \frac{x^{2n+1}}{(2n+1)!} = \sum_{n = 0}^{\infty} B_{2n+1}\left(\frac{1}{3}\right) \frac{x^{2n+1}}{(2n+1)!}
 \end{align*}
which implies,
\begin{align*}
 B_{2n+1}\left(\frac{1}{3}\right) =  -\frac{(2n+1)I_{2n}}{3^{2n+1}}.
 \end{align*}
