"""
This script illustrates that, in a cyclic cubic field, the index of the group of cyclcotomic units in the group of units is a multiple of the class number.
"""

"""
Let K be a number field. This function takes as parameters h an element in K (and implicitely deduces K), epsilons which is a base of the unit group quotient by the torsion subgroup, rs which is a list of complex roots of a defining polynomial of K. If h is a unit, the function computes exponents e1,...,er such that

 h = epsilon1^e1 * epsilon2^e2 * ... * epsilonr^er.

"""
def exponents(h,epsilons,rs):
    M=Matrix([[log(epsj.polynomial()(ri).abs()) for epsj in epsilons] for ri in rs])
    b=Matrix([[log((h.polynomial()(ri)).abs())] for ri in rs ]) 
    sol=M.solve_right(b)
    l=list(sol.transpose()[0])
    return map(round,l)

"""
Given a matrix M, compute its echelon form Me and compute a basis of a vector space which complement the image of M.
"""
def img_and_complement(M):
    rs=M.echelon_form().rows()
    n=M.ncols()
    zero = n*[0]
    Me=Matrix([ro for ro in rs if list(ro) != zero])
    k=Me.nrows()
    ii=[i for i in range(k) if Me[i][i] == 0] + range(k,n)
    return Me,ii 

"""
Given a in (Z/mZ)^* return the corresponding automorphism of the m-th cyclotomic field.
"""
def morph(a):
    def sigma(z):
        Qm=z.parent()
        zm=Qm.gen()
        return z.polynomial()(zm^a)
    return sigma

"""
Given an odd integer m, the function returns the list qs of prime powers dividing m and a generator for (Z/q^eZ)* for each prime power q^e in qs.
"""
def gens_mod(m,ell=3):
    mfact=m.factor()
    gs=[]
    qs=[]
    for qe in m.factor():
       q,e=qe
       qi=q^e
       qs.append(qi)
       gi=IntegerModRing(q^e).multiplicative_generator().lift()
       gi_=crt([gi,1],[q^e,m//(q^e)])
       gs.append(gi_) 
    return gs,qs

"""
Let m be the m-th cyclotomic field Qm and f a polynomial defining a subfield K of cyclic Galois group of degree p. Compute a set of representatives for Gal(K/Q) where Q is the field of rationals.
"""
def representatives(f,Qm):
    thetaf=f.roots(Qm)[0][0]
    zm=Qm.gen()
    K.<aK>=NumberField(f)
    m=K.disc().sqrt()
    ell = f.degree()
    gs,qs=gens_mod(m)
    n=len(gs)
    ggs=[]
    for i in range(n):
        gi=IntegerModRing(qs[i])(gs[i])
        oi=gi.multiplicative_order()
        ggi=(gi^(oi//ell)).lift()
        ggs.append(ggi)
    import itertools
    rowsM=[]
    for a_ in itertools.product(range(ell),repeat=n):
        a = crt([gs[i]^list(a_)[i] for i in range(n)],qs)
        if morph(a)(thetaf) == thetaf:
            rowsM.append(list(a_))
    Me,ii=img_and_complement(Matrix(GF(ell),rowsM))
    hs=[gs[i] for i in ii]
    cartesian=[]
    for i in range(len(gs)):
        ggi = ggs[i]
        oo = multiplicative_order(IntegerModRing(m)(gs[i]))
        cartesian.append([ggi^j for j in range(oo//ell)])
    reps = []
    for g_ in itertools.product(*cartesian):
        reps.append(crt(list(g_),qs) % m)
    for r in Me.rows():
        tmp=[]
        a=prod([gs[i]*r[i].lift() for i in range(n)])
        a= a%m
        for i in range(ell):
            tmp = tmp + [rep*a^i % m for rep in reps]
        reps = tmp[:]
    return hs,reps 
     
"""
Given a polynomial f which defines a cyclic cubic field, compute the exponents of the cyclotomic units with respect to a system of fundamental units. 
"""
def exponents_from_poly(f):
    K.<a_>=NumberField(f)
    OK=K.ring_of_integers()
    m=ZZ(K.disc().sqrt())
    Qm.<zm>=NumberField(cyclotomic_polynomial(m))
    ell = f.degree()
    m1=euler_phi(m) // ell
    hs,reps = representatives(f,Qm)
    g1=hs[0] # only for cyclic fields
    sigma=morph(g1) 
    u2=IntegerModRing(m)(1/2).lift()
    eta0=(zm^(g1*u2)-zm^(-g1*u2))/(zm^u2-zm^(-u2))
    assert eta0.is_unit()
    eta=prod([morph(a)(eta0) for a in reps])
    assert eta != 1
    thetaf=f.roots(Qm)[0][0]
    K.<theta>=NumberField(f)
    eps0_=K.units()[0]
    eps1_=K.units()[1]
    eps0=eps0_.polynomial()(thetaf)
    eps1=eps1_.polynomial()(thetaf)
    epss=[eps0,eps1]
    zmRR=cyclotomic_polynomial(m).complex_roots()[0]
    ell = f.degree()
    rs=[zmRR^(g1^i % m) for i in range(0,ell-1)]
    exps=exponents(Qm(eta),epss,rs)
    sigma_m_eps0=sigma(eps0)
    Mat=Matrix(ZZ,[exponents(eta,epss,rs),exponents(sigma(eta),epss,rs)])
    return Mat

Qx.<x>=QQ['x']
f=x^3-7*x^2-10*x-1
print exponents_from_poly(f)