/*
 * Implementation of the formula for the number of subgroups of an abelian group
 * given the paper G.A Miller "On the Subgroups of an Abelian Group",
 * The Annals of Mathematics, 2nd Ser., Vol. 6, No. 1 (Oct., 1904), pp. 1-6
 * http://dx.doi.org/10.2307/2007151
 *
 * ver. 1.1 (c) 2007 by Max Alekseyev
 */


\\ number of subgroups of the group C_{p^a1} x C_{p^a2} x ... C_{p^ak} where 
\\ p is prime and a=[a1,a2,...,ak] with 1 <= a1 <= a2 <= ... <= ak
{ numsubgrp(p,a,b=0) = local(r=0,q,n,m,n1,m1);
  a = vecsort(a);
  if(b==0,
    forvec(c=vector(#a,i,[0,a[i]]), 
      if(c==0,next); 
      r+=numsubgrp(p,a,c), 
    1);
    return(r + 1);
  );

  \\ removing leading zeros from b
  q=1; while(q<=#b && b[q]==0,q++);
  if(q>#b,return(0));
  b = vecextract(b,2^#b-2^(q-1));

  \\ initialize q, m[], n[]
  q = #b;
  m = n = vector(b[q]);
  m1 = 1;
  n1 = 1; 
  for(i=1,b[q],
    while(i>a[m1],m1++);
    m[i] = #a - m1 + 1;
    while(i>b[n1],n1++);
    n[i] = q - n1 + 1;
  );

  prod(i=1,q, (p^(m[b[q-i+1]]-i+1) - 1) / (p^(n[b[q-i+1]]-i+1) - 1) ) * 
  p^sum(x=0,q-1, (q-x) * sum(j=if(x<1,1,b[x]),b[x+1]-1, m[j] - n[j]) )
}