"""
Computing the Alexander polynomial of a SnapPy manifold inside Sage.

Example:

----------------------------------------------------------------------
| Sage Version 4.0, Release Date: 2009-05-29                         |
| Type notebook() for the GUI, and license() for information.        |
----------------------------------------------------------------------
sage: import snappy
sage: from alexander import alexander_polynomial
sage: M = Manifold('m004')
sage: alexander_polynomial(M)
-a^2 + 3*a - 1
sage: M = snappy.Manifold('v3551')
sage: alexander_polynomial(M)
a^4*b^4 + a^4*b^3 + a^3*b^4 + a^4*b^2 + 2*a^3*b^3 + a^2*b^4 + 2*a^3*b^2 + 2*a^2*b^3 + a^3*b + 2*a^2*b^2 + a*b^3 + 2*a^2*b + 2*a*b^2 + a^2 + 2*a*b + b^2 + a + b + 1
sage: M = snappy.Manifold('9^4_1')
sage: alexander_polynomial(M)
a^2*c*d - a^2*b*d - 2*a*c*b*d + 2*a*b^2*d - 2*a^2*c + a^2*b + 2*a*c*b - a*b^2 - a*c*d + 2*a*b*d + c*b*d - 2*b^2*d + 2*a*c - 2*a*b - c*b + b^2

By Nathan M. Dunfield.   This code is released into the public domain.  
"""

import os, sys, re, string
import snappy
import sage
from sage.all import Integers, vector, matrix, LaurentPolynomialRing, PolynomialRing, gcd, prod


#----------------------------------------------------------------
#
#  Abelianization of the fundamental group
#
#----------------------------------------------------------------

ZZ = Integers()

def abelianize_word(word, gens):
    return vector(ZZ, [ word.count(g) - word.count(g.swapcase()) for g in gens])

class MapToFreeAbelianization(sage.structure.sage_object.SageObject):
    def __init__(self, fund_group):
        self.domain_gens = fund_group.generators()
        R = matrix(ZZ, [abelianize_word(R, self.domain_gens) for R in fund_group.relators()]).transpose()
        D, U, V = R.smith_form()
        self.U = U
        self.elementry_divisors = [D[i,i] for i in range(D.ncols())] + [0,]*(D.nrows() - D.ncols())

    def range(self):
        return ZZ**self.elementry_divisors.count(0)

    def __call__(self, word):
        D = self.elementry_divisors
        v = self.U*abelianize_word(word, self.domain_gens)
        return vector(ZZ, [v[i] for i in range(len(D)) if D[i] == 0])


#--------------------------------------------------------------
#
#  Computing the Alexander polynomial
#
#--------------------------------------------------------------

class MapToGroupRingOfFreeAbelianization(MapToFreeAbelianization):
    def __init__(self, fund_group):
        MapToFreeAbelianization.__init__(self, fund_group)
        n = self.elementry_divisors.count(0)
        self.R = LaurentPolynomialRing(ZZ, list(string.ascii_lowercase[:n]))

    def range(self):
        return self.R

    def __call__(self, word):
        v = MapToFreeAbelianization.__call__(self, word)
        return prod([ g**v[i] for i, g in enumerate(self.R.gens())])

def fox_derivative(word, phi, var):
    R = phi.range()
    if len(word) == 0:
        return R(0)
    w = word[0]
    if w == w.lower():
        term = R(0) if w != var else R(1)
    else:
        term =   R(0) if w.lower() != var else -phi(var.upper())
    return term + phi(w) * fox_derivative(word[1:], phi, var)

# It's best to deal with matrixes of polynomials rather than Laurent
# polynomials, so we need to be able to clear denominators.  This add
# to the complexity of the code.  

def join_lists(list_of_lists):
    ans = []
    for L in list_of_lists:
        ans += L
    return ans

def uniform_poly_exponents(poly):
    return [list(e) if hasattr(e, "__getitem__") else (e,) for e in poly.exponents()]
    
def minimum_exponents(elts):
    A =  matrix(ZZ, join_lists([ uniform_poly_exponents(p) for p in elts]))
    return vector( [min(row) for row in A.transpose()] )

def convert_laurent_to_poly(elt, expshift, P):
   if elt == 0:
       return P(0)
   return sum( [  c*prod([g**e for g, e in zip(P.gens(), vector(exps) + expshift)]) for c, exps in zip(elt.coefficients(), uniform_poly_exponents(elt))])

def alexander_polynomial(manifold):
    G = manifold.fundamental_group()
    phi = MapToGroupRingOfFreeAbelianization(G)
    R = phi.range()
    P = PolynomialRing(R.base_ring(), R.gens_dict().keys())
    M = [[fox_derivative(rel, phi, var)  for rel in G.relators()] for  var in G.generators()]
    expshift = -minimum_exponents(join_lists(M))
    M = matrix(P, [[ convert_laurent_to_poly(p, expshift, P) for p in row] for row in M])
    alex_poly = gcd(M.minors(G.num_generators() - 1))
    # Normalize it
    return convert_laurent_to_poly(alex_poly, -minimum_exponents( [alex_poly] ), P)