Question

Let's say you want to compute the Hall polynomial $g^\lambda_{(r),\mu}(p)$. According to [Dutta and Prasad, Degenerations and orbits in finite abelian groups], the orbits under the automorphism group of a finite abelian group are given by $\{O_I|I\subset J(P_\lambda)\}$, where $J(P_\lambda)$ denotes the lattice of order ideals in a certain poset $P_\lambda$. The same paper also gives a formula for $|O_I|$.

Clearly, the Hall polynomial that you are looking for is $(p^r-p^{r-1})^{-1}\sum_{I} |O_I|$, the sum being over all $I$ for which the order of an element in $O_I$ is $p^r$ and for which the quotient of the group of type $\lambda$ by any element of $O_I$ is of type $\mu$. As pointed out by David Speyer in the comments to an earlier version of this answer, $I$ is uniquely determined by these conditions. So a final form of the answer is obtained by explaining how to obtain $I$ from $\lambda$ and $\mu$.

Given an element $(p^{v_1},p^{v_2},\dotsc)$, the type of the quotient is found by computing the Smith canonical form of the matrix $\begin{pmatrix} p^{\lambda_1} & & &\\ & p^{\lambda_2} & &\\ & & \ddots & \\& & &p^{\lambda_l}\\ p^{v_1} & p^{v_2} & \cdots & p^{v_l} \end{pmatrix}$.

By the characterization of order ideals in $P_\lambda$, we have that $v_i\leq v_{i-1}\leq v_i+(\lambda_{i-1}-\lambda_i)$.

Proposition. Let $I\subset P_\lambda$ be an order ideal. Let $\mu$ be the type of the group obtained by going modulo an element of $O_I$. Then

$\mu_l=v_l$

$\mu_{l-1}=\lambda_l+v_{l-1}-v_l$

$\mu_{l-2}=\lambda_{l-1}+v_{l-2}-v_{n-1}$

$\vdots$

$\mu_{i}=\lambda_{i+1}+v_i-v_{i+1}$

$\vdots$

$\mu_1=\lambda_2+v_1-v_2$.

Proof. The gcd of $i\times i$ minors of the above matrix can be seen to be $v_{l-i+1}+\lambda_{l-i}+\dotsb+\lambda_l$ v_{l-i+1}+\lambda_{l-i+2}+\dotsb+\lambda_l$ (using the inequalities on $v_i$). Therefore, we get $\mu_{l-i+1}+\dotsb+\mu_l=v_{l-i+1}+\lambda_{l-i}+\dotsb+\lambda_l$, \mu_{l-i+1}+\dotsb+\mu_l=v_{l-i+1}+\lambda_{l-i+2}+\dotsb+\lambda_l$, from which the above identities follow.QED.

This allows us to recover $(v_1,v_2,\dotsc)$ once we know $\lambda$ and $\mu$, where (this also follows from the above proposition), $\mu$ is a partition obtained from $\lambda$ by removing a horizontal strip of length $r$. In particular, $v$ is uniquely determined by $\lambda$ and $\mu$.

We get

$v_l=\mu_l$, and $v_i=\mu_i-[(\lambda_{i+1}+\dotsb+\lambda_l)-(\mu_{i+1}+\dotsb+\mu_l)]$ for $i<l$.

Let $I\subset P_\lambda$ be the order ideal defined by $(v_1,v_2,\dotsc)$. We get

$g^\lambda_{(r)\mu}(p)=|O_I|/(p^r-p^{r-1})$.

It easily follows (from the formula for $|O_I|$ in our paper, which says that $|O_I|$ is a monic polynomial in $p$ of degree $\sum_i (\lambda_i-\nu_i)$) that $g^\lambda_{(r)\mu}(p)$ is monic in $p$ of degree $n(\lambda)-n(\mu)$. This could perhaps give another approach to Hall's theorem (in analogy with the proof that MacDonald gave).