p,q=4,1

F1.<a1,a2,a3,a4,a5,a6,a7,a8,f> = FreeGroup()
S.<t> = LaurentPolynomialRing(ZZ)
G1 = F1 / [a1/a6/a1*a3, a3/a8/a3*a1, a2/a5,a4/a7,     # Wirtinger relation
           a1*a2*a3*a4*f^-4,                          # lens relation
           a1^-1*f*a1/f*a1/a5,
           (a1*a2)^-1*f*a1*a2/a1/f*a1*a2/a6,
           (a1*a2*a3*a4)^-1*f*a1*a2*a3*(a1*a2)^-1/f^-1*a1*a2*a3*a4/a7,
           (a1*a2*a3*a4)^-1*f*a1*a2*a3*a4*(a1*a2*a3)^-1/f*a1*a2*a3*a4/a8]
alex_matrix=Ge.alexander_matrix([t,t,t,t,t,t,t,t,1])
alex_poly=gcd(ae.minors(8))
alex_matrix
print "Alexander poly:", alex_poly

OUTPUT:

[  -t^-1 + 1           0        t^-1           0           0          -1           0           0           0]
[       t^-1           0   -t^-1 + 1           0           0           0           0          -1           0]
[          0           1           0           0          -1           0           0           0           0]
[          0           0           0           1           0           0          -1           0           0]
[          1           t         t^2         t^3           0           0           0           0       4*t^4]
[          1           0           0           0          -1           0           0           0    t^-1 - 1]
[          0           1           0           0           0          -1           0           0 t^-2 - t^-1]
[          0           0        t^-1   -t^-1 + 1           0           0          -1           0 t^-4 + t^-3]
[          0           0           0           1           0           0           0          -1 t^-4 - t^-3]
Alexander poly: t^-6