Formula or code to compute number of subgroups of a certain order of an abelian $p$-group

Given a finite abelian $p$-group and its factorization into groups of the form $\mathbb{Z}/p^k\mathbb{Z}$, does anyone know of a formula that gives the number of subgroups of a certain index/order? As I'm sure such a formula would contain some nasty product or sum, is there a computer algebra system out there that knows how to compute this?

-
I don't know if this will help, but [this question/answer] gives a way of describing all subgroups of a direct sum of cyclic groups of prime power order, in terms of certain congruences. – Arturo Magidin Dec 7 '11 at 5:22
I probably need to code this up myself, since I've needed this a lot in the stuff that I've been working on lately. Somehow I think that needing this can't be too rare, so I would expect some CAS to have implemented this. :) – pki Dec 7 '11 at 6:38
Check the GAP documentation, or ask in the GAP mailing list. – Arturo Magidin Dec 7 '11 at 6:40
In answer to a similar question on MO, mathoverflow.net/questions/78956/…, Greg Martin gave a formula for the number of subgroups of a specific isomorphism type, so you could sum over the possible isomorphism types to get what you want. – Derek Holt Dec 7 '11 at 9:52

Let $\alpha = (\alpha_1,\dots,\alpha_\ell)$ be a partition, so that $\alpha_1\ge\cdots\ge\alpha_\ell$. (In this formula it is convenient to allow some of the parts of the partition at the end to equal 0.) Define the notation $${\mathbb Z}_\alpha = {\mathbb Z}/p^{\alpha_1}{\mathbb Z} \times \cdots \times {\mathbb Z}/p^{\alpha_\ell}{\mathbb Z}$$ for a general $p$-group of type $\alpha$. Define similarly a partition $\beta$, and suppose that $\beta\preceq\alpha$, meaning that $\beta_j\le\alpha_j$ for each $j$. We want to count the number of subgroups of ${\mathbb Z}_\alpha$ that are isomorphic to ${\mathbb Z}_\beta$.
Let $a=(a_1,\dots,a_{\alpha_1})$ be the conjugate partition to $\alpha$, so that $a_1=\ell$ for example; similarly, let $b$ be the conjugate partition to $\beta$. Then the number of subgroups of ${\mathbb Z}_\alpha$ that are isomorphic to ${\mathbb Z}_\beta$ is $$\prod_{i=1}^{\alpha_1} \genfrac{[}{]}{0pt}{}{a_i-b_{i+1}}{b_i-b_{i+1}}p^{(a_i-b_i)b_{i+1}},$$ where $$\genfrac{[}{]}{0pt}{}nm = \prod_{j=1}^m \frac{p^{n-m+j}-1}{p^j-1}$$ is the Gaussian binomial coefficient.
To answer your specific question, you'd want to sum over subpartitions $\beta\preceq\alpha$ such that $\beta_1$ equals the exponent in question.