In the paper under review, the authors consider the following generating function \[ V_t(q)=\sum_{n\geq 0} M(t,n)q^n:=\sum_{1\leq k_1\leq k_2\leq \cdots \leq k_t}\frac{q^{k_1+k_2+\cdots+k_t}}{(1-q^{k_1})^2 (1-q^{k_2})^2\cdots (1-q^{k_t})^2}. \] Note that when \(t = 1\), this again is the generating function \(\sum_{n\geq 0} \sigma_1(n)q^n\) for the sum of the divisors of \(n\). The authors reveal a rich arithmetic aspect of these generalized divisor functions. For instance, it was proven that \[ M(2, 5n + 1) \equiv 0\pmod 5 \] and many other similar identities in Theorem 8.1. In addition to those surprising congruences, the authors also obtained a novel identity among these generalized divisor sums. Namely, \[ V_t(q)=\sum_{1\leq k_1\leq k_{2}\leq \cdots \leq k_{2t-1}}\frac{k_1q^{k_1+k_2+\cdots+k_{2t-1}}}{(1-q^{k_1}) (1-q^{k_2})\cdots (1-q^{k_{2t-1}})}. \] The main approach to prove the above identity involves rational function approximation to MacMahon-type generating functions. One such example involves multiple \(q\)-harmonic sums \[ \sum\limits_{k=1}^n\frac{(-1)^{k-1}\begin{bmatrix} n \\ k\end{bmatrix}_q(1+q^k)q^{\binom{k}{2}+tk}}{[k]_q^{2t}\begin{bmatrix} n+k\\ k\end{bmatrix}_q} =\sum\limits_{1\leq k_1\leq\cdots\leq k_{2t}\leq n}\frac{q^{n+k_1+k_3\cdots +k_{2t-1}}+q^{k_2+k_4+\cdots +k_{2t}}}{[n+k_1]_q[k_2]_q\cdots [k_{2t}]_q}. \]