#This program has been checked with $\\texttt{SageMath 10.5}$ with $\\texttt{Sirocco}$ so that the fundamental group package can run. Otherwise, a presentation for the groups has been included for convenience. The following code computes twisted Alexander polynomials.
##################################################################
#Twisted Fox derivative of x_n with respect to x_m
#associated with the morphism "hom" (hom(x_i)=hom[i]) and
#the linear representation "subs_repres".
#OPTIONS: "mmod" can be either a number or a list [n,p1,...]
# where n is a natural number and p1,... is a list of
# one-variable polynomials
def tw_dif(n,subs_repres,hom,m,R):
    QR=R.fraction_field()
    QR.inject_variables(verbose=False)
    subs_repres_aux=[_.change_ring(QR) for _ in subs_repres]
    dd=subs_repres_aux[0].dimensions()
    if dd[0]!=dd[1]:
        return 'Wrong representation matrix'
    else:
        d=dd[0]
    uno=identity_matrix(d)
    cero=0*uno
    mod_rels=[]
    if m<=0:
        return 'Error'
    elif n==-m:
        return (-subs_repres_aux[m-1]^(-1))*(1/hom[m-1])
    elif n==m:
        return uno
    else:
        return cero
#################################################################
#Twisted Fox derivative of "list" with respect to x_m
#associated with the morphism "hom" (hom(x_i)=hom[i]) and
#relative to the linear representation "subs_repres"
#OPTIONS: "mmod" can be either a number or a list [n,p1,...]
# where n is a natural number and p1,... is a list of
# one-variable polynomials
def tw_dif_list(list,subs_repres,hom,m,R):
    QR=R.fraction_field()
    QR.inject_variables(verbose=False)
    l=len(list)
    subs_repres_aux=[_.change_ring(QR) for _ in subs_repres]
    if l==1:
        return tw_dif(list[0],subs_repres_aux,hom,m,R)
    elif sign(list[0])==-1:
        aux=subs_repres_aux[abs(list[0])-1]^(-1)*(1/hom[abs(list[0])-1])
    else:
        aux=subs_repres_aux[abs(list[0])-1]*hom[abs(list[0])-1]
    return tw_dif(list[0],subs_repres_aux,hom,m,R)+aux*tw_dif_list([list[i+1] for i in range(l-1)],subs_repres_aux,hom,m,R)

########################################################################
#tw_alex_matrix(list_rels,subs_repres,hom) it produces the Alexander
#matrix for the twisted case relative to the relations "list_rels",
#to the morphism "hom" and to the linear representation "subs_repres".
#OPTIONS: "mmod" can be either a number or a list [n,p1,...]
# where n is a natural number and p1,... is a list of
# one-variable polynomials
def tw_alex_matrix(list_rels,subs_repres,hom,R):
    QR=R.fraction_field()
    QR.inject_variables(verbose=False)
    subs_repres_aux=[_.change_ring(QR) for _ in subs_repres]
    dd=subs_repres_aux[0].dimensions()
    if dd[0]!=dd[1]:
        return 'Wrong representation matrix'
    else:
        d=dd[0]
    lrels=len(list_rels)
    lrepres=len(subs_repres_aux)
    m=max(flatten(list_rels))
    if m!=lrepres:
        return print("error")
    A_aux=block_matrix(lrels,m,[tw_dif_list(i,subs_repres_aux,hom,j,R) for i in list_rels for j in [1..m]])
    A_aux_simpl=[QR(_) for _ in A_aux.list()]
    A=matrix(dd[0]*lrels,dd[0]*m,A_aux_simpl)
    return A

########################################################################
#rep_check(As,rels) checks if the list of matrices "As" produces a representation
#of a finitely presented group whose relations are given in "rels".
def rep_check(As,rels):
    n=len(As)
    d=As[0].dimensions()[0]
    for k in As:
        if k.determinant()!=1:
            print(k)
            return 'False'
    for k in rels:
        for i in range(n):
            kk=identity_matrix(d)
            for j in range(len(k)):
                kk=kk*As[abs(k[j])-1]^(sign(k[j]))
            if kk!=identity_matrix(d):
                print(j,kk)
                return 'False'
    return 'True'

########################################################################
#reps(As,rels) produces a list of equations required by the entries of the matrices "As"
#to define a representation of a finitely presented group whose relations are given in "rels".
def reps(As,rels):
    n=len(As)
    d=As[0].dimensions()[0]
    eqs_aux=[_.det()-1 for _ in As]
    for k in rels:
        for i in range(n):
            kk=identity_matrix(d)
            for j in range(len(k)):
                kk=kk*As[abs(k[j])-1]^(sign(k[j]))
            kk=kk-identity_matrix(d)
            eqs_aux=eqs_aux+[_.numerator() for _ in kk.list()]
    return eqs_aux

#Define the equations of the curves $C_1$ and $C_2$ and calculate the fundamental groups of their complements. If
#$\texttt{Sirocco}$ is not installed, skip this cell.
var('x,y,z')
A.<x,y> = (QQ[I])[]
f=x^6+y^6+6*x*y
g1aux=(f+4)*(f-4)
g2aux=(f+4)*(f-4*I)
g1=A(g1aux)
g2=A(g2aux)
C1=Curve(g1)
C2=Curve(g2)
G1=C1.fundamental_group()
G2=C2.fundamental_group()
G1_rels=[list(_.Tietze()) for _ in G1.relations()]
G2_rels=[list(_.Tietze()) for _ in G2.relations()]

G1_rels=[[-3, -2, -1, 2, 3, -2, 1, 2], [1, 2, -3, -2, 3, -1, 2, 3, -2, -2, -1, 2], [-1, -3, -1, 3, -2, -1, 2, 1, 3, 2, 3, -2], [3, 1, 2, 2, 3, -2, -1, -3, -2, -2, -1, 2], [2, 3, 2, -3, -2, -2, 1, 2, -3, 1, -2, -1], [-1, 2, 3, 3, -2, -3, -2, -2, 3, -2, -1, 2, 2, 1, 2, -1]]
G2_rels=[[1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2]]"

#We check that both Characteristic Varieties coincide
#G1.characteristic_varieties()
#G2.characteristic_varieties()

#Define the number field $F=\\mathbb Q[\\xi]$, where $\\xi$ is an 8-primitive root of unity and $R=F[t,t^{-1}]$
#Laurent polynomial ring over $F$.

var('t')
sqI = -QQ[x](x^4+1).roots(QQbar)[0][0]
F = NumberField((sqI).minpoly(), 'sqI', embedding=sqI)
F.inject_variables(verbose=False)
R= LaurentPolynomialRing(F,[t], order='degrevlex')
R.inject_variables(verbose=False)

#Define a $U(2)$ representation for $G_1$.
A1=matrix(2,2,[sqI,0,0,-sqI^3]);
A2=matrix(2,2,[-(sqI^3 - sqI)/2,-sqI^2/2,-sqI^2,-(sqI^3 - sqI)/2]);
A3=matrix(2,2,[-(sqI^3 - sqI)/2,-1/2,1,-(sqI^3 - sqI)/2]);
As=[A1,A2,A3];

#Check As is a valid representation, that is, the matrices of As satisfy the relations in the presentation of $G_1$.
rep_check(As,G1_rels)

#Calculate the twisted Alexander polynomial associated with the representation As anf the morphism
#$\varepsilon:G_1\to \mathbb Z$ defined as $\varepsilon(x_1)=\varepsilon(x_2)=\varepsilon(x_3)=1$.
A=tw_alex_matrix(G1_rels,As,[t,t,t],R)
A_red=A.delete_columns([0,1])
i=4
Fi=[_.numerator() for _ in A_red.minors(i)];
Delta1=gcd(Fi)
Delta0=(As[0]*t-identity_matrix(2)).determinant()
(Delta1/Delta0).factor()