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?

As pointed out in the comments, this question was answered by Greg Martin on MathOverflow. I include their answer below for completeness:

I had to look this up as well at some point in my research. The answer is yes, and a Google search for "number of subgroups of an abelian group" leads to several downloadable papers, not all of them easy to read. The paper "On computing the number of subgroups of a finite abelian group" by T. Stehling, in Combinatorica 12 (1992), contains the following formula and (I think) references to where it has appeared earlier in the literature.
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.

