Can this nested sum be expressed in terms of generalized harmonic numbers and the cycle index polynomials of the symmetric groups? For a paper I was working on recently I needed to find the value of the following sum:
$$S(n,k) = \sum_{i_1 = 1}^n \sum_{i_2 = i_1+1}^n \cdots \sum_{i_k=i_{k-1}+1}^n \frac{1}{i_1 i_2 \cdots i_k}.$$ 
I found a couple of references (by Adamchik and Cheon and El-Mikkaway) that have an expression for $S(n,k)$ as a polynomial containing generalized harmonic numbers $H_n^{(r)}$, where
$$H_n^{(r)} = \sum_{j=1}^n \frac{1}{j^r}.$$
For example, 
$$S(n,2) = \frac{1}{2}\left(H_n^2 - H^{(2)}_n \right),$$
$$S(n,3) = \frac{1}{6}\left(H_n^3 - 3H_n H^{(2)}_n + 2 H_n^{(3)}\right),$$
$$S(n,4) = \frac{1}{24}\left(H_n^4 - 6 H^2_n H_n^{(2)} + 3 (H_n^{(2)})^2 + 8 H_n H_n^{(3)}  - 6 H_n^{(4)}\right).$$
Neither of these papers considers the corresponding polynomial sequence of indeterminates (the polynomials before substituting in the generalized harmonic numbers), though.  Calculations for small values of $n$ indicate that these are the cycle index polynomials of the symmetric groups, with the sign pattern such that each factor of $H_n^{(r)}$ contributes a $+1$ if $r$ is odd and a $-1$ if $r$ is even.  
Could someone give a proof of this, particularly one with a combinatorial flavor that gives some real insight into why the cycle index polynomials of the symmetric groups show up here (assuming that they do)?   
(For the record, I don't need this answered for my paper.  I just want to know for my own sake.)
As a side note, the papers also give the extremely (and, to me, surprisingly) simple expression $$S(n,k) = \frac{1}{n!} \left[ n+1 \atop k+1 \right],$$ where $\left[ n \atop k \right]$ is an unsigned Stirling number of the first kind.
 A: I'll edit in the details later, but this is essentially a consequence of Polya's enumeration theorem together with an inclusion-exclusion argument; see this blog post.
Edit:  Okay, so it's a little easier than that; this is just a well-known identity for symmetric functions
$$e_n(x_1, x_2, ...) = \frac{1}{n!} \sum_{\sigma \in S_n} \text{sgn}(\sigma) p_{\sigma}$$
in disguise, where $p_{\sigma} = p_{\lambda_1} ... p_{\lambda_k}$ if $\sigma$ has cycle type $(\lambda_1, ... \lambda_k)$ and $p_k = \sum_i x_i^k$.  You get your identity upon substituting $x_i = \frac{1}{i}$.  The identity above is equivalent to the generating function identity
$$\prod (1 + x_i t) = \exp \left( p_1 t - \frac{p_2 t^2}{2} + \frac{p_3 t^3}{3} \mp ... \right)$$
which formally follows from the identity
$$\prod \frac{1}{1 - x_i t} = \exp \left( p_1 t + \frac{p_2 t^2}{2} + \frac{p_3 t^3}{3} + ... \right)$$
which is in turn a consequence of Polya's enumeration theorem.  These are more refined versions of the identities which appear in a previous MO answer of mine, and both of the arguments in that answer prove this identity without change (you just have to keep track of all the weights); it's a weighted sum over injective functions.  (These identities also automatically imply the Stirling number identity.)
A: One can rephrase Qiaochu's argument without using symmetric functions as follows:
A standard property of $\left[ n+1\atop k+1\right]$ is the generating function
  $$ \sum_{k=0}^{n+1}\left[ n+1\atop k\right]t^k = t(t+1)(t+2)\cdots (t+n). $$
Equating coefficients of $t^{k+1}$ gives the result
   $$ S(n,k) =\frac{1}{n!}\left[ n+1\atop k+1\right]. $$
Now if the cycle index of the symmetric group $S_k$ is $Z_k(x_1,x_2,\dots,x_k)$,
then it is well-known that
  $$ \sum_k Z_k(x_1,x_2,\dots)t^k  = \exp\sum_{r\geq 1}x_r\frac{t^r}{r}. $$
Put $x_r=(-1)^{r-1}H_n^{(r)}$ to get
  $$ \sum_{k\geq 0} Z_k(H_n^{(1)},-H_n^{(2)},\dots)t^k = \left(1+t\right)\left(1+\frac t2\right)\cdots         
     \left(1+\frac tn\right), $$
from which it is immediate that $S(n,k)=Z_k(H_n,-H_n^{(2)},\dots,(-1)^{k-1}H_n^{(k)})$.
A: Are the formulae for $S(n,m)$ above in terms of the generalised Harmonic numbers $H_n(r)=\sum_{k=1}^n\frac{1}{k^r}$ actually correct?
Mathematica, for example, shows up $S(10,2)=1026576$ (for Stirling number of the first kind), whereas $\frac{1}{2}\left(H_{10}^2-H_{10}^{(2)}\right) = \frac{177133}{50400}$.
Further, assuming the Pochhammer symbol is defined $(a)_n=a(a+1)(a+2)\cdots(a+n-1)$ then Adamchik's recursive formula is incorrect too!
I'm actually trying to get to grips with these formulae myself, so it's quite frustrating I can't get them to work. Perhaps I'm using incorrect conventions for symbols? However, I have checked and double checked everything. Can anyone help?
A: Correction: the first paragraph is nonsense; I got confused with the notation for $S(n,m)$ for Stirling numbers, and the finite sum form (it was a late night...)
However, if you follow Adamchik's formula from his original paper you do get
$S(n,2) = 1/2(H_n - H_n^{(2)})$,
i.e. the power of $H_n$ is equal to $1$, and not equal to $2$, as expected. So the second paragraph still holds from what I can see. Any help on this matter would be much appreciated.
