1 files changed, 320 insertions, 252 deletions

diff --git a/kummer_degree.sage b/kummer_degree.sage
index 161e8c9..3085efd 100644
--- a/kummer_degree.sage
+++ b/kummer_degree.sage
@@ -1,16 +1,3 @@
-#############################################################################
-# This program allows one to compute the degree of certain field extensions #
-# of the rational numbers. In particular, it can compute the degree over Q  #
-# of extensions of the form Q( sqrt[N](G), \zeta_M ), where:                #
-#  - N and M are integers with N dividing M;                                #
-#  - G is a finitely generated subgroup of the multiplicative group of Q    #
-#  - \zeta_M is a primitive M-th root of unity                              #
-# The group G does not have to be given in a particular format. A finite    #
-# set of generators is sufficient.                                          #
-#############################################################################
-
-from sage.misc.cachefunc import CachedFunction
-
 # Computes the "adelic Kummer failure", i.e. the degrees of the intersection
 # of the the Kummer extension Q(\sqrt{2^n}{G}) with the M-th cyclotomic field
 # over Q_{2^n}.
@@ -18,12 +5,12 @@ from sage.misc.cachefunc import CachedFunction
 # Input: a good basis B for the torsion-free group G, organized as a list of
 # lists, and a non negative integer d. They have to satisfy the following:
 # 1. Each list B[i] contains all basis elements of 2-divisibility i.
-# 2. The basis given by B is 2-maximal, that is to say it satisfies Theorem 
-#    14 of Debry-Perucca; i.e. each element of B[i] is, up to plus or minus 
-#    1, the 2^i-th power of a strongly 2-indivisible rational.
+# 2. The basis given by B is 2-maximal, that is to say it satisfies Theorem 14
+#    of Debry-Perucca; i.e. each element of B[i] is, up to plus or minus 1, the
+#    2^i-th power of a strongly 2-indivisible rational.
 # 3. For i != d every element of B[i] is positive, and B[d][0] is the only
 #    negative element of B[d] (if d=-1, then there is no negative element).
-# 3'. Notice that the existence on negative elements in B[0] does not change
+# 3'. Notice that the existence on negative elements in B[0] does not influence
 #    the correctness of the algorithm; in fact, the function adjust_sign
 #    produces a basis that may have negative elements of divisibility 0, and
 #    this basis is given as input for adelic_failure_gb.
@@ -32,9 +19,9 @@ from sage.misc.cachefunc import CachedFunction
 # integer such that the intersection of Q(\sqrt{2^n}{G}) with Q_\infty is
 # contained in Q_M0.
 # The table ad_fail has N rows, where N is defined below.
-# Each row R=ad_fail[i] contains a variable number of pairs (d,r), where d 
-# is a divisor of M0 and r is the degree of Q(\sqrt{2^{i+1}}{G}) \cap Q_d 
-# over Q_{2^n}.
+# Each row R=ad_fail[i] contains a variable number of pairs (d,r), where d is a
+# divisor of M0 and r is the degree of Q(\sqrt{2^{i+1}}{G}) \cap Q_d over
+# Q_{2^n}.
 # Each divisor of M appears at most once on each row, and the last element of
 # the last row is of the form (M0,r0).
 def adelic_failure_gb( B, d ):
@@ -42,48 +29,48 @@ def adelic_failure_gb( B, d ):
     # The table to be returned (or printed at the end), as described above.
     ad_fail = []
 
-    # N is such that for every n > N the adelic failure of Q(\sqrt{2^n}{G})
-    # is the same as that of Q(\sqrt{2^N}{G}).
-    # We always have to include n=3, because of problem with sqrt(2) in Q_8
-    # (in theory, this is not necessary in some cases, e.g. if 2 does not 
-    # divide any element of G).
-    # If the negative generator is on the last level, we need to increase N
-    # by 1, because it would contribute to the shortlist in the next level
-    # (by taking the root of an even power).
+    # N is such that for every n > N the adelic failure of Q(\sqrt{2^n}{G}) is
+    # the same as that of Q(\sqrt{2^N}{G}).
+    # We always have to include n=3, because of problem with sqrt(2) in Q_8 (in
+    # theory, this is not necessary in some cases, e.g. if 2 does not divide
+    # any element of G).
+    # If the negative generator is on the last level, we need to increase N by
+    # 1, because it would contribute to the shortlist in the next level (by
+    # taking the root of an even power).
     if d == len(B)-1:
         N = max(3,len(B)+1)
     else:
         N = max(3,len(B))
 
     # The intersection is given by adding the square roots of the elements of
-    # this shortlist (and a "special element", not always of the form sqrt d,
-    # coming from taking a suitable root of a negative generator; this
-    # special element is dealt with later). The shortlist grows at each step,
-    # so we declare it before starting to loop over n and we build it
-    # incrementally. The "special element" is of the form \zeta_{2^n}\sqrt b,
-    # which we encode as (n,b). We use the value (1,1) (special element -1)
-    # to say that there is no special element at this level.
+    # this shortlist (and a "special element", not always of the form sqrt(d),
+    # coming from taking a suitable root of a negative generator; this special
+    # element is dealt with later). The shortlist grows at each step, so we
+    # declare it before starting to loop over n and we build it incrementally.
+    # The "special element" is of the form \zeta_{2^n}\sqrt{b}, which we encode
+    # as (n,b). We use the value (1,1) (special element -1) to say that there
+    # isno special element at this level.
     shortlist = []
     special_element = (1,1)
 
     # The integers M, giving the smallest cyclotomic field in which lies the
-    # whole intersection with Q_\infty, grows with n. As with the shortlist,
-    # we declare it here and increase it appropriately at each step.
+    # whole intersection with Q_\infty, also grows with n. As with the
+    # shortlist, we declare it here and increase it appropriately at each step.
     M = 1
 
-    for n in range( 1, N+1 ):
+    for n in range( 1, N+1 ): # 1 \leq n \leq N
 
         # We add the new elements to the shortlist, modifying M if needed.
-        # This is not done in case we are in the extra "fake" level (this
-        # case dealt with immediately below).
+        # This is not done in case we are in the extra "fake" level (this case
+        # dealt with immediately below).
         if n-1 < len(B):
             for g in B[n-1]:
                 # Case of negative g
                 if g < 0 and n > 1:
                     # Special element of the form \zeta_{2^{n+1}}\sqrt(b).
-                    # It is contained in Q_{lcm(2^{n+1},cyc_emb(b))}, except 
-                    # in case n=2 and b=2, in which case it is contained in
-                    # Q_4. We store it as as pair ( n+1, b ).
+                    # It is contained in Q_{lcm(2^{n+1},cyc_emb(b))}, except in 
+                    # case n=2 and b=2, in which case it is contained in Q_4.
+                    # We store it as as pair ( n+1, b ).
                     special_element = ( n+1, abs(g)^(1/(2^(n-1))) )
                     M = lcm( M, special_embed( special_element ) )
                 else:
@@ -91,8 +78,8 @@ def adelic_failure_gb( B, d ):
                     shortlist.append( b )
                     M = lcm( M, cyc_embed(b) )
 
-        # We add a root of an even power of the negative generator, as soon
-        # as we are beyond its level.
+        # We add a root of an even power of the negative generator, as soon as
+        # we are beyond its level.
         if d != -1 and n == d+2:        
             b = abs(B[d][0])^(1/2^d)
             shortlist.append( b ) 
@@ -114,20 +101,19 @@ def adelic_failure_gb( B, d ):
 
         # For each divisor dM of M, compute the degree of the intersection of
         # Q(\sqrt{2^n}{G}) with Q_dM over Q_{2^n}. We need to compute the
-        # number r of elements in the subgroup of G generated by the
-        # shortlist that lie in Q_dM. This is going to be a power of 2. The
-        # sought degree will be r, up to considering some special cases
-        # (described below).
+        # number r of elements in the subgroup of G generated by the shortlist
+        # that lie in Q_dM. This is going to be a power of 2. The sought degree
+        # will be r, up to considering some special cases (described below).
         # 
-        # This algorithm could be inefficient for groups of big rank. It can
-        # be improved by precomputing subgroups of G/G^2 and the cyclotomic
-        # fields containing them.
+        # This algorithm could be inefficient for groups of big rank. It can be
+        # improved by precomputing subgroups of G/G^2 and the cyclotomic fields
+        # containing them.
 
         aux = [] # Next line of ad_fail table
 
         for dM in divisors( M ):
-            # We only care of the intersection with Q_dM if it contains the
-            # 2^n roots of unity.
+            # We only care of the intersection with Q_dM if it contains the 2^n
+            # roots of unity.
             if dM % (2^n) != 0:
                 continue
 
@@ -144,7 +130,7 @@ def adelic_failure_gb( B, d ):
             if n <= d and dM % (2^(n+1)) == 0 and n > 1:
                 r *= 2
 
-            # We loose a factor of 2 if we have sqrt(2) in Q_8 or \zeta_8 in
+            # We loose a factor of 2 if we have sqrt(2) in Q_8 or \zeta_8 2 in
             # Q_4.
             if 8 in H and dM % 8 == 0 and (n >= 3 or (n == 2 and n <= d)):
                 r = r/2
@@ -153,11 +139,11 @@ def adelic_failure_gb( B, d ):
             # We have to consider all possible special elements arising from
             # multiplying the given element with the other elements of the
             # shortlist.
-            # If any of them is of the form \zeta_8 2q for q a square, there
-            # is nothing to do: in fact \zeta_8 2 embeds in Q_4, which
-            # coincides with Q_{2^n} (n must be 2); so we would double the 
-            # degree because of the existence of the special element, but 
-            # then we would loose another factor of 2 because of this.
+            # If any of them is of the form \zeta_8 2q for q a square, there is
+            # nothing to do: in fact \zeta_8 2 embeds in Q_4, which coincides
+            # with Q_{2^n} (n must be 2); so we would double the degree because
+            # of the existence of the special element, butthen we would loose
+            # another factor of 2 because of this.
             if special_element != (1,1) and special_element[0] == n+1:
                 nothing_to_do = False
                 intersecting_QdM = False
@@ -184,8 +170,11 @@ def cyc_embed( b ):
         m *= 4
     return abs(m)
 
-# Computes the minimal cyclotomic field containing \zeta_{2^n}\sqrt(b).
-def special_embed( (n,b) ):
+# Computes the minimal cyclotomic field containing \zeta_{2^n}\sqrt(b),
+# where (n,b)=p
+def special_embed( p ):
+    n = p.first
+    n = p.second
     m = squarefree_part(b)
     if n == 3 and m % 2 == 0:
         return 4 * cyc_embed(m/2)
@@ -201,7 +190,7 @@ def special_embed( (n,b) ):
 #     degree of Q_{l^m,l^n} over Q_{l^m}, where l = L[i]
 def total_l_adic_failure( B ):
 
-    bp = bad_primes(B)
+    bp = bad_primes( B )
     ret = ( bp, [] )
 
     for l in bp:
@@ -210,32 +199,35 @@ def total_l_adic_failure( B ):
     return ret
 
 # Computes the "bad primes", i.e. the ones for which the l-adic part is not
-# maximal. Used by total_l_adic_failure.
+# maximal. Used by l_adic_degree and total_l_adic_failure.
 def bad_primes( B ):
     M = exponent_matrix( B )
     (a,b) = M.dimensions()
-
     if a > b or M.rank() < a:
-        print "This is not a basis"
+        print( "This is not a basis" )
         return
 
     # Compute which primes l divide all minors of the exponent matrix
-    bad_p = list( prime_factors( gcd( M.minors(a) ) ) )
-    if 2 not in bad_p:
-        bad_p += [2] # 2 is always bad
+    ms = M.minors( a )
+    d = ms[0]
+    for m in ms:
+        d = gcd( d, m )
+    bad_primes = list( prime_factors( d ) )
+    if 2 not in bad_primes:
+        bad_primes += [2] # 2 is always bad
+    bad_primes.sort() # Ensures 2 is always first
 
-    return sorted(bad_p) # Ensures 2 is always first
+    return bad_primes
    
-# Computes the l-adic failure for a specific l. Returns a "table" as
-# described above "total_l_adic_failure".
+# Computes the l-adic failure for a specific l. Returns a "table" as described
+# above "total_l_adic_failure".
 # B is any basis for G.
 def l_adic_failure( B, l ):
 
     r = len(B)
     GB = make_good_basis( B, l )
     
-    # Computes the parameters for G over Q4. For l odd, they are the same as
-    # over Q.
+    # Computes the parameters over Q4. For l odd, they are the same as over Q.
     p = parameters_Q4( GB, l )
     maxM = max( [ sum(x) for x in p ] )
     maxN = max( maxM, len(GB)-1 )
@@ -277,6 +269,7 @@ def l_adic_failure_from_data( B, l, tablel, M, N ):
         return 1
     r = len(B)
 
+    # Basically, the "dual" of what we do in l_adic_degree.
     if m > len(tablel):
         if n > len(tablel):
             return l^tablel[-1][1]
@@ -285,70 +278,127 @@ def l_adic_failure_from_data( B, l, tablel, M, N ):
     
     return l^(r*n-tablel[m-1][0][n-1])
 
-# Returns the l-dvisibility parameters of G over Q, given a good basis b of
-# G, as a list of pairs (di,hi).
-def parameters_Q( b, l ):
-    if l != 2:
-        return [x for a in [[(i,0)]*len(b[i]) for i in range(len(b))] \
-                for x in a]
-    else:
-        return [x for a in [[(i,1-max(0,sgn(g))) for g in b[i]] \
-                for i in range(len(b))] for x in a]
-
-# Computes the l-divisibility parameters of G over Q4, given a good basis gb
+# Computes the l-divisibility parameters of G over Q4, given a good basis b
 # over Q for G. Returns a list of pairs (di,hi).
 # If l is odd it just uses the good basis given to compute the parameters.
-# If l=2, it uses the results of PST-2.
 def parameters_Q4( gb, l ):
+    # Converts from "good basis format" to simple list
+    b = []
+    for x in gb:
+        b += x
+    ret = []
+    
     if l != 2:
-        return parameters_Q( gb, l )
+        for i in range( len( gb ) ):
+            for j in gb[i]:
+                ret.append( (i,0) )
+        return ret
     else:
-        # Converts from "good basis format" to simple list
-        b = [ x for a in gb for x in a ]
-    
-        M = exponent_matrix( b )
+        R.<y> = PolynomialRing( QQ )
+        pol = R(y^2+1)
+        Q4.<eye> = NumberField( pol ) # I already use i for other things
 
-        # Thanks to my bad notation, the elements of b are what are called
-        # g_i in the article, while the elements of bb will be the b_i's.
-        bb = [ abs(b[i])^(2^(-divisibility(M[i],2))) for i in range(len(b)) ]
+        # Factorize basis elements over Q4 and so on.
+        d = []
+        B = []
+        h = []
+        ideals_list = set()
+        M = [] # Exponent matrix of the Bi's
 
-        # Computing a combination of bb elements of the form 2 * square.
-        MM = exponent_matrix( bb + [2] ).change_ring( GF( 2 ) ) 
-        
-        for a in MM.kernel().basis():
-            if a[-1] != 0:
-                # the vector a[0:-1] gives the coefficients for a combination
-                # of the b_i's of the form 2 * square.
-                
-                # Index of asis elements that do appear in the combination
-                c = [ i for i in range(len(a[0:-1])) if a[i] != 0 ]
-                
-                # Change of basis to include the element of the form 
-                # 2 * square to some power.
-                div_max = -1
-                i_max = -1
-                for j in c:
-                    if divisibility(M[j],2) > div_max:
-                        div_max = divisibility(M[j],2)
-                        i_max = j
-                b[i_max] = prod([b[j]^(2^(div_max-divisibility(M[j],2))) \
-                                 for j in c])
-                M = exponent_matrix(b)
+        # Pre-process to find all ideals appearing in the factorization and fix
+        # a chosen generator for each of them. This is important in order to
+        # compute the "sign" (h-parameter) of an element with respect to it Bi.
+        for g in b:
+            factorization_list = list( Q4.ideal(g).factor() )
+            ideals_list |= set( [ x[0] for x in factorization_list ] )
+        ideals_list = list( ideals_list )
+        # Chooses a generator of each principal ideal in the list
+        irreducibles_list = [ J.gens_reduced()[0] for J in ideals_list ]
+
+        # Compute the Q4-parameters of the given basis b. Also computes the
+        # exponent matrix of the Bi's 
+        for g in b:
+            factorization_list = list( Q4.ideal(g).factor() )
+            exps = [ x[1] for x in factorization_list ]
+            d.append( divisibility( exps, l ) )
+            Bg = 1
+            for j in range(len(factorization_list)):
+                a = 0
+                for i in range( len( ideals_list ) ):
+                    if ideals_list[i] == factorization_list[j][0]:
+                        a = irreducibles_list[i]
+                        break
+                Bg *= a ^ (exps[j]/(l^d[-1]))
+            B.append(Bg)
+            u = g / (Bg^(l^d[-1]))
+            if not u.is_unit():
+                print( "Error: g is not the right power of the computed Bg." )
+                print( "g:", g, ", Bg:", Bg, ", exponent:", l^d[-1] )
+            if u == 1:
+                h.append( 0 )
+            elif u == -1:
+                h.append( 1 )
+            else:
+                h.append( 2 ) 
+
+        # Make the exponent matrix M (for now as a list of rows)
+        for g in B:
+            row = [0] * len(ideals_list)
+            for i in range(len(ideals_list)):
+                I = ideals_list[i]
+                ee = 1
+                while (I^ee).divides(g):
+                    ee += 1
+                row[i] = ee-1
+            M.append(row)
 
-                # Return the parameters, changing only those of the element
-                # of highest divisibility that appears in the combination.
-                ret = [(divisibility(M[i],2),1-max(0,sgn(b[i]))) \
-                        for i in range(len(b)) ]
-                d1, h1 = ret[i_max]
-                if d1 == 1:
-                    h1 = 1 - h1
-                elif d1 == 0:
-                    h1 = 2
-                ret[i_max] = ( d1+1, h1 )
+        # If the Bi's are not strongly independent, apply the algorithm (only
+        # once) to produce a new basis. The new basis has maximal parameters.
+        coeffs = find_combination( matrix(M), l )
+        if coeffs != []:
 
-                return ret
+            maxi = -1
+            maxd = -1
+            for i in range(len(d)):
+                if d[i] > maxd and coeffs[i] != 0:
+                    maxd = d[i]
+                    maxi = i
 
-        return parameters_Q( gb, 2 )
+            x = [(a/coeffs[maxi]).lift() for a in coeffs] # Now a vector of int
+
+            new_element = 1
+            for i in range(len(d)):
+                new_element *= b[i]^( x[i] * l^(d[maxi]-d[i]) )
+        
+            b[maxi] = new_element
+
+            # Compute new B, d and so on.
+            factorization_list = list( Q4.ideal(b[maxi]).factor() )
+            exps = [ x[1] for x in factorization_list ]
+            d[maxi] = divisibility( exps, l )
+            Bg = 1
+
+            for j in range(len(factorization_list)):
+                a = 0
+                for i in range( len( ideals_list ) ):
+                    if ideals_list[i] == factorization_list[j][0]:
+                        a = irreducibles_list[i]
+                        break
+                Bg *= a ^ (exps[j]/(l^d[maxi]))
+            B[maxi] = Bg
+            M[maxi] = [ x[1] for x in list( Q4.ideal(Bg).factor() ) ]
+            u = b[maxi] / (Bg^(l^d[maxi]))
+            if not u.is_unit():
+                print( "Error: new element is not the right power of B." )
+                print( "New el.:",b[maxi],", B:",Bg,", exponent:",l^d[maxi] )
+            if u == 1:
+                h[maxi] = 0
+            elif u == -1:
+                h[maxi] = 1
+            else:
+                h[maxi] = 2 
+            
+        return [(d[i],h[i]) for i in range(len(d))]
 
 # Uses Theorem 18 to compute the degree of Kummer extensions.
 def compute_vl( p, n, m, r ):
@@ -358,26 +408,16 @@ def compute_vl( p, n, m, r ):
     
     return M - m + r*n - sum( ni )
 
-# Computes a basis for G from a frozen_set of generators.
-# Returns a pair (B,torsion) where B is a list of rational numbers and
-# torsion is true if -1 is in G and false otherwise.
-@CachedFunction
-def generators_to_basis( G ):
-    (BM,BM_primes) = exponent_matrix_with_sign_and_primes( list(G) )
-    BM = BM.echelon_form()
-    BM_primes.append(-1)
-
-    B = [prod([BM_primes[i]^r[i] for i in range(len(r))]) for r in BM.rows()]
-
-    return ( [ x for x in B if abs(x) != 1 ], -1 in B )
-
-
 # Given any basis b of a group G computes an l-good basis for G. This is done
 # using the algorithm outlined in the proof of Theorem 14 of (Debry-Perucca).
 def make_good_basis( b, l ):
     M = exponent_matrix( b )
-    d = [ divisibility(x,l) for x in M ]
-    B = [ abs(b[i])^(1/(l^d[i])) for i in range(len(b)) ]
+    d = []
+    B = []
+    for i in range(len(b)):
+        di = divisibility( M[i], l )
+        d.append( di )
+        B.append( abs(b[i])^(1/(l^di)) )
 
     # Computes the coeffiecients of a linear combination of the rows of M
     # that is zero modulo l. These coefficients are elements of F_l.
@@ -385,8 +425,8 @@ def make_good_basis( b, l ):
 
     while coeffs != []:
 
-        # Computes which basis element (with non-zero coefficient in the
-        # linear combination above) has maximal divisibility.
+        # Computes which basis element (with non-zero coefficient in the linear
+        # combination above) has maximal divisibility.
         maxi = -1
         maxd = -1
         for i in range(len(d)):
@@ -394,24 +434,31 @@ def make_good_basis( b, l ):
                 maxd = d[i]
                 maxi = i
 
-        x = [ (a/coeffs[maxi]).lift() for a in coeffs ] # Now a vector of int
+        x = [ (a/coeffs[maxi]).lift() for a in coeffs ] # Now a vector of ints
 
-        new_elt = prod([b[i]^(x[i]*l^(d[maxi]-d[i])) for i in range(len(d))])
+        new_element = 1
+        for i in range(len(d)):
+            new_element *= b[i]^( x[i] * l^(d[maxi]-d[i]) )
         
-        b[maxi] = new_elt
+        b[maxi] = new_element
         M = exponent_matrix( b )
         d[maxi] = divisibility( M[maxi], l )
         B[maxi] = abs(b[maxi])^(1/(l^d[maxi]))
         
         coeffs = find_combination( exponent_matrix( B ), l )
 
-    return [[b[i] for i in range(len(b)) if d[i]==j] for j in range(max(d)+1)]
+    GB = [[]]
+    for i in range(len(d)):
+        while( len(GB) <= d[i] ):
+            GB.append([])
+        GB[d[i]].append(b[i])
+    return GB
 
-# Takes a good basis B and adjusts the sign of the elements so that there is
-# at most one negative generator (of positive divisibility). The input is a
-# good basis in the format returned by make_good_basis.
-# Returns a pair (B,d), where B is the updated basis and d is the 
-# divisibility parameter of the only negative element remained.
+# Takes a good basis B and adjusts the sign of the elements so that there is at
+# most one negative generator (of positive divisibility). The input is a good
+# basis in the format returned by make_good_basis.
+# Returns a pair (B,d), where B is the updated basis and d is the divisibility
+# parameter of the only negative element remained.
 # The sign of the d=0 elements is just ignored in the other steps of the
 # algorithm, so we keep them negative.
 def adjust_sign( B ):
@@ -430,7 +477,6 @@ def adjust_sign( B ):
 # Given the exponent matrix M of a list of rational numbers B, returns the
 # coefficients of a linear combination that is weakly l-divisible, or [] if
 # the B[i] are strongly l-independent.
-@CachedFunction
 def find_combination( M, l ):
     M = M.change_ring( GF( l ) )
     if M.rank() != min( M.dimensions() ):
@@ -443,8 +489,8 @@ def find_combination( M, l ):
 # Returns the minimal l-valuation of the exponents.
 def divisibility( A, l ):
     if len(A) == 0:
-        print "Warning: computing the divisibility of a torsion element.",
-        print "Returning +Infinity."
+        print( "Warning: computing the divisibility of a torsion element.", end="" )
+        print( "Returning +Infinity." )
         return +Infinity
     return min( [ valuation( x, l ) for x in A ] )
 
@@ -469,10 +515,11 @@ def exponent_matrix( B ):
 # Returns a pair (M,L) where M is the modified exponent matrix and L is the
 # list of primes appearing in the factorization.
 def exponent_matrix_with_sign_and_primes( B ): 
-
-    prime_list = list({x for a in [prime_factors(g) for g in B] for x in a})
+    prime_list = set()
+    for g in B:
+        prime_list |= set( prime_factors( g ) )
+    prime_list = list( prime_list )
     np = len( prime_list )
-
     rows = []
     for g in B:
         rowg = [0] * np
@@ -480,19 +527,21 @@ def exponent_matrix_with_sign_and_primes( B ):
             for i in range( np ):
                 if f[0] == prime_list[i]:
                     rowg[i] = f[1]
-        rowg.append( 1 - max(0,sgn(g)) ) 
+        s = 0
+        if sgn(g) == -1:
+            s = 1
+        rowg.append( s )
         rows.append( rowg )
     return ( matrix( rows ), prime_list )
 
-# The function TotalKummerFailure (with its many wrapper functions) computes
-# the total table of Kummer Failures for the given group G. Moreover, the
-# result is cached for speeding up following computations on the same group,
-# even if the same group is given with a different set of generators.
-#
+# This is a wrapper function for total_kummer_failure( G, True ), see below.
+def TotalKummerFailure( G ):
+    total_kummer_failure( G, True )
+
 # Input: any set of generators for a subgroup G of Q*.
 # If output=False, returns a 4-uple (t,MM,NN,D):
-# - t is a pair, where t[0] is the rank of G and t[1] is either True (if G
-#   has torsion) or False (is it does not).
+# - t is a pair, where t[0] is the rank of G and t[1] is either True (if G has
+#   torsion) or False (is it does not).
 # - MM is the pair (M0,divisors(M0))
 # - NN is the pair (N0,divisors(N0))
 # - D is a table F, where F[j][i] is the ration between phi(m)n^r and the
@@ -501,34 +550,42 @@ def exponent_matrix_with_sign_and_primes( B ):
 #   computed for a torsion-free part of G.
 # If output is True, outputs this data in a human-readable way and does not
 # return any value.
+def total_kummer_failure( G, output ):
 
-def TotalKummerFailure( G ):
-    total_kummer_failure( G, True )
-
-@CachedFunction
-def total_kummer_failure_cachable( G ):
-    B, torsion = generators_to_basis( G )
-    return total_kummer_failure_cachable_basis( tuple(B), torsion )
+    # Computing a basis.
+    (BM,BM_primes) = exponent_matrix_with_sign_and_primes( G )
+    BM = BM.echelon_form()
+    BM_primes.append(-1)
+    B = []
+    torsion = False
 
-@CachedFunction
-def total_kummer_failure_cachable_basis( B, torsion ):
-    B = list(B)
-    r = len(B) # Rank of G
+    for r in BM.rows():
+        gr = product( [ BM_primes[i]^r[i] for i in range(len(r)) ] )
+        if gr == -1:
+            torsion = True
+            break
+        elif gr == 1:
+            break
+        else:
+            B.append(gr)
 
-    if r == 0:
-        print "G is torsion. The extension is cyclotomic. Stopping."
+    if len(B) == 0:
+        print( "G is torsion. The extension is cyclotomic. Stopping." )
         return False
 
+    r = len(B) # Rank of G
+
     # Compute l-adic data (straightforward)
-    ( bad_p, l_adic_failure_table ) = total_l_adic_failure( B )
+    ( bad_primes, l_adic_failure_table ) = total_l_adic_failure( B )
 
     # Compute adelic data.
     (GB,d) = adjust_sign( make_good_basis( B, 2 ) )
     adelic_failure_table = adelic_failure_gb( GB, d )
     
     # Computing the bounds M0 and N0
-    N0 = prod( [ bad_p[i] ^ len( l_adic_failure_table[i][-1][0] ) \
-            for i in range(len(bad_p)) ] )
+    N0 = 1
+    for i in range(len( bad_primes )):
+        N0 *= bad_primes[i] ^ len( l_adic_failure_table[i][-1][0] )
     # Extra factors of 2 may come from the adelic failure
     N0 = lcm( N0, 2^len( adelic_failure_table ) )
     divs_N0 = divisors(N0)
@@ -555,20 +612,16 @@ def total_kummer_failure_cachable_basis( B, torsion ):
                                 FT[j][l] = lcm( FT[j][l], pp[1] )
 
     # Adding l-adic failure to the table
-    for i in range( len( bad_p ) ):
-        l = bad_p[i]
+    for i in range( len( bad_primes ) ):
+        l = bad_primes[i]
         for j in range(len(divs_N0)):
             dN = divs_N0[j]
             fl = l_adic_failure_from_data(B,l,l_adic_failure_table[i],dN,dN)
             for h in range(len(divs_M0)):
                 FT[j][h] *= fl
 
-    return ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT )
+    ret = ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT )
 
-def total_kummer_failure( G, output ):
-    
-    ret = total_kummer_failure_cachable( frozenset(G) )
-    
     if output:
         print_total_table( ret )
         # Uncomment following line for case list description.
@@ -581,9 +634,9 @@ def total_kummer_failure( G, output ):
 # as the ration between 2^eN^r and the degree of Q_{M,N} over Q_M, where e=1
 # if N is even and e=0 otherwise.
 
-# Makes the failure table for the torsion case when M/N is even. In this
-# case an entry of the table is doubled if the corresponding value of N
-# (actually, of gcd(N,N0) ) is even, and is kept the same otherwise.
+# Makes the failure table for the torsion case when M/N is even. In this case,
+# an entry of the table is doubled if the corresponding value of N (actually,
+# of gcd(N,N0) ) is even, and is kept the same otherwise.
 # The expected degree (over Q) 2^e * phi(M) * N^r, where e=1 if N is even and
 # e=0 otherwise.
 def torsion_table_even( data ):
@@ -600,8 +653,8 @@ def torsion_table_even( data ):
 
 # Makes the failure table for the torsion case when M/N is odd. In this case
 # the entry at (M,N) is taken from the torsion-free entry at (2M,N).
-# In other words, the expected degree (over Q) 2^e * phi(M) * N^r, where e=1
-# if N is even and e=0 otherwise.
+# In other words, the expected degree (over Q) 2^e * phi(M) * N^r, where e=1 if
+# N is even and e=0 otherwise.
 def torsion_table_odd( data ):
 
     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
@@ -625,101 +678,114 @@ def print_total_table( data ):
 
     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
     
-    print "M_0 =", M0
-    print "N_0 =", N0
-    print ""
-    print "The following table shows the total failure of Kummer", 
-    if torsion:
-        print "degrees in case the quotient M/N is EVEN."
-    else:
-        print "degrees."
-    print "The degree of the Kummer extension (M,N) is e / f, where",
+    print( "M_0 =", M0 )
+    print( "N_0 =", N0 )
+    print( "" )
+    print( "The following table shows the total failure of Kummer degrees", end="" )
     if torsion:
-        print "e = phi(M)*N^rank(G) if N is odd and e = 2*phi(M)*N^rank(G)",
-        print "if N is even",
+        print( "in\n case the quotient M/N is EVEN." )
     else:
-        print "e = phi(M)*N^rank(G)",
-    print "and f is the entry of the table below at the row labelled with",
-    print "gcd(N,N0) and the column labelled with gcd(M,M0)."
-    print ""
-
+        print( "." )
+    print( "Columns correspond to values of M, rows to values of N" )
+    print( "" )
+    print( "The degree of the Kummer extension (M,N) can be extracted by taking" )
+    print( "the value f (failure) of the entry at (gcd(N,N0),gcd(M,M0)) and" )
+    print( "simply computing ed(M,N) / f, where ed(M,N) is the expected degree" )
+    print( "of the Kummer extension." )
     if torsion:
+        print( "In this case (-1 is in G), we have ed(M,N) = 2^e*phi(M)*N^r," )
+        print( "where e=1 if N is even and e=0 if N is odd." )
         FT1 = torsion_table_even( data )
     else:
+        print( "In this case (G is torsion-free) we have ed(M,N) = phi(M)*N^r," )
         FT1 = FT
+    print( "where r is the rank of G." )
+    print( "" )
 
     tt = [ ["","|"] + divs_M0 ]
     tt.append( "-" * (len(divs_M0)+2) )
     for i in range(len(divs_N0)):
         tt.append( [ divs_N0[i] ] + ["|"]  + FT1[i] )
-    print table(tt)
-    print ""
+    print( table(tt) )
+    print( "" )
 
     if torsion:
-        print "The following table shows the total failure of Kummer degrees",
-        print "if the quotient M/N is ODD and is read as the previous one."
-        print ""
+        print( "The following table shows the total failure of Kummer degrees in" )
+        print( "case the quotient M/N is ODD." )
+        print( "This table can be read exactly as the first one." )
+        print( "" )
 
         # A good strategy is the following:
         #   A little translation exercise: the failure at (M,N) in the torsion
-        #   case is either the same as that for (2M,N) in the torsion-free
-        #   case (if M is even) or its double (if M is odd).
+        #   case is either the same as that for (2M,N) in the torsion-free case
+        #   (if M is even) or its double (if M is odd).
         # However, due to problems in reading the table for bigger M, it is
         # easier to just compute the degree every time, and then deduce the
-        # failure. This is not too inefficient, since we can use the data
-        # that we have already computed via kummer_degree_from_total_table.
+        # failure. This is not too inefficient, since we can use the data that
+        # we have already computed via kummer_degree_from_total_table.
         new_FT = torsion_table_odd( data )
         # Printing the new table
         tt = [ ["","|"] + divs_M0 ]
         tt.append( "-" * (len(divs_M0)+2) )
         for i in range(len(divs_N0)):
             tt.append( [ divs_N0[i] ] + ["|"]  + new_FT[i] )
-        print table(tt)
-        print ""
+        print( table(tt) )
+        print( "" )
 
 def print_case_list( data ):
 
     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
-    FT_odd = torsion_table_odd( data )
 
+    FT_odd = torsion_table_odd( data )
+    # FT1 is either FT or FT_even in the torsion case
+    FT1 = FT
     if torsion:
         FT1 = torsion_table_even( data )
-        pf = sorted(list({ x for row in FT_odd for x in row }))
-    else:
-        FT1 = FT
-        pf = sorted(list({ x for row in FT1 for x in row }))
-
+    pf = []
+    for row in FT1:
+        pf += row
+    if torsion:
+        for row in FT_odd:
+            pf += row
+    pf = list(set(pf))
+    pf.sort()
     for f in pf:
-        print "Failure is", f, "if",
+        print( "Failure is", f, "if", end="" )
         if torsion:
-            print "M/N is EVEN and",
-        print "(gcd(M,M0),gcd(N,N0)) is one of the following:"
-        print [ (divs_M0[j], divs_N0[i]) \
-                for i in range(len(divs_N0)) for j in range(len(divs_M0)) \
-                if FT1[i][j] == f ]
-
+            print( "M/N is EVEN and", end="" )
+        print( "(gcd(M,M0),gcd(N,N0)) is one of the following:" )
+        lijst = []
+        for i in range(len(divs_N0)):
+            for j in range(len(divs_M0)):
+                if FT1[i][j] == f:
+                    lijst.append( ( divs_M0[j], divs_N0[i] ) )
+        print( lijst )
         if torsion:
-            print "or if M/N is ODD and (gcd(M,M0)),gcd(N,N0)) is one of the",
-            print "following:"
-            print [ (divs_M0[j], divs_N0[i]) \
-                    for i in range(len(divs_N0)) for j in range(len(divs_M0))\
-                    if FT_odd[i][j] == f ]
+            print( "or if M/N is ODD and (gcd(M,M0)),gcd(N,N0)) is one of the", end="" )
+            print( "following:" )
 
-        print ""
+            lijst_odd = []
+            for i in range(len(divs_N0)):
+                for j in range(len(divs_M0)):
+                    if FT_odd[i][j] == f:
+                        lijst_odd.append( ( divs_M0[j], divs_N0[i] ) )
+            print( lijst_odd )
+
+        print( "" )
 
 # Extracts a specific value of failure from the total table.
 def kummer_failure_from_total_table( M, N, data ):
     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
-
+    FT1 = FT
     if torsion:
         if (M/N) % 2 == 0:
             FT1 = torsion_table_even( data )
         else:
             FT1 = torsion_table_odd( data )
-    else:
-        FT1 = FT
 
-    return FT1[divs_N0.index( gcd( N, N0 ) )][divs_M0.index( gcd( M, M0 ) )]
+    i = divs_N0.index( gcd( N, N0 ) )
+    j = divs_M0.index( gcd( M, M0 ) )
+    return FT1[i][j]
 
 # Computes the degree of the Kummer extension (M,N), by taking as input the
 # table computed by TotalKummerFailure.
@@ -737,17 +803,20 @@ def kummer_degree_from_total_table( M, N, data ):
 # M must be a multiple of N.
 def KummerDegree( G, M, N ):
     if M % N != 0:
-        print "M is not a multiple of N"
+        print( "M is not a multiple of N" )
         return -1
 
     data = total_kummer_failure(G,False)
     ((r,torsion),(M0,divs_M0),(N0,divs_N0),FT) = data
 
-    e_deg = euler_phi(M) * N^r
+    exp_deg = euler_phi(M) * N^r
     
     if torsion and N % 2 == 0:
-        e_deg *= 2
+        exp_deg *= 2
     
+    j = divs_M0.index(gcd(M,M0))
+    i = divs_N0.index(gcd(N,N0))
+
     if torsion: 
         if (M/N)%2 == 0:
             failure = torsion_table_even( data )
@@ -756,5 +825,4 @@ def KummerDegree( G, M, N ):
     else:
         failure = FT
 
-    return e_deg/failure[divs_N0.index(gcd(N,N0))][divs_M0.index(gcd(M,M0))]
-
+    return exp_deg / failure[i][j]