## USAGE:
## Start sage in command line. Then:
## sage: load("Abelian_lattice_count.py")
## sage: LC([1,50])
## 93
## sage: LC([4,4,4])
## 11287

## INPUT: LC([d_1, ... ,d_n])
## OUTPUT: number of polarized isogenies
## from a (d_1, ... , d_n) polarized abelian variety
## to a principally polarized abelian variety.
## The output number equals \nu(d_1, ..., d_n) defined in Section 1.3 of
## "Curve counting on abelian varieties" (with J. Bryan, R. Pandharipande, and Q. Yin)

## For information on individual subgroups, use LC([d_1, .. , d_n], True):
## sage: LC( [1,6] , True)
## Count for G = Multiplicative Abelian group isomorphic to C6
## Subgroup H = Multiplicative Abelian subgroup isomorphic to C2 x C3 generated by {f} ; Order of gens [6] ; Contribution 6
## Subgroup H = Multiplicative Abelian subgroup isomorphic to C3 generated by {f^2} ; Order of gens [3] ; Contribution 3
## Subgroup H = Multiplicative Abelian subgroup isomorphic to C2 generated by {f^3} ; Order of gens [2] ; Contribution 2
## Subgroup H = Trivial Abelian subgroup ; Order of gens [] ; Contribution 1
## 12

## By Georg Oberdieck, June 2015.

def LC( d , verbose = False):
  #d = [d_1, d_2, ... , d_n]
  assert(x > 0 for x in d)
  d = [x for x in d if x != 1]

  G = AbelianGroup(len(d), d)
  if verbose: print "Count for G =", G
  res = 0
  for H in G.subgroups():
    h = H.gens()
    a = [x.order() for x in h]
    r = len(h)
    res += prod(  number_divisors( a[i], a[j] )  for (i,j) in xmrange([r,r]) if i >= j )
    if verbose: print "Subgroup H =", H, "; Order of gens", a, "; Contribution", prod(  number_divisors( a[i], a[j] )   for (i,j) in xmrange([r,r]) if i >= j )
  return res

@cached_function
def number_divisors(a, b):
  #returns number of elements x in Z/b s.t. ax = 0 (b).
  #This equals #Hom(Z/a, Z/b)
  return sum( 1 for m in range(b) if m*a % b == 0)

##END