kummer-degrees

Compute the degree of Kummer extensions
git clone https://git.tronto.net/kummer-degrees
Download | Log | Files | Refs | README | LICENSE

kummer_degree.sage (30746B)


      1 # Computes the "adelic Kummer failure", i.e. the degrees of the intersection
      2 # of the the Kummer extension Q(\sqrt{2^n}{G}) with the M-th cyclotomic field
      3 # over Q_{2^n}.
      4 #
      5 # Input: a good basis B for the torsion-free group G, organized as a list of
      6 # lists, and a non negative integer d. They have to satisfy the following:
      7 # 1. Each list B[i] contains all basis elements of 2-divisibility i.
      8 # 2. The basis given by B is 2-maximal, that is to say it satisfies Theorem 14
      9 #    of Debry-Perucca; i.e. each element of B[i] is, up to plus or minus 1, the
     10 #    2^i-th power of a strongly 2-indivisible rational.
     11 # 3. For i != d every element of B[i] is positive, and B[d][0] is the only
     12 #    negative element of B[d] (if d=-1, then there is no negative element).
     13 # 3'. Notice that the existence on negative elements in B[0] does not influence
     14 #    the correctness of the algorithm; in fact, the function adjust_sign
     15 #    produces a basis that may have negative elements of divisibility 0, and
     16 #    this basis is given as input for adelic_failure_gb.
     17 #
     18 # Output: a table ad_fail as described below. Call M0 the smallest positive
     19 # integer such that the intersection of Q(\sqrt{2^n}{G}) with Q_\infty is
     20 # contained in Q_M0.
     21 # The table ad_fail has N rows, where N is defined below.
     22 # Each row R=ad_fail[i] contains a variable number of pairs (d,r), where d is a
     23 # divisor of M0 and r is the degree of Q(\sqrt{2^{i+1}}{G}) \cap Q_d over
     24 # Q_{2^n}.
     25 # Each divisor of M appears at most once on each row, and the last element of
     26 # the last row is of the form (M0,r0).
     27 def adelic_failure_gb( B, d ):
     28 
     29     # The table to be returned (or printed at the end), as described above.
     30     ad_fail = []
     31 
     32     # N is such that for every n > N the adelic failure of Q(\sqrt{2^n}{G}) is
     33     # the same as that of Q(\sqrt{2^N}{G}).
     34     # We always have to include n=3, because of problem with sqrt(2) in Q_8 (in
     35     # theory, this is not necessary in some cases, e.g. if 2 does not divide
     36     # any element of G).
     37     # If the negative generator is on the last level, we need to increase N by
     38     # 1, because it would contribute to the shortlist in the next level (by
     39     # taking the root of an even power).
     40     if d == len(B)-1:
     41         N = max(3,len(B)+1)
     42     else:
     43         N = max(3,len(B))
     44 
     45     # The intersection is given by adding the square roots of the elements of
     46     # this shortlist (and a "special element", not always of the form sqrt(d),
     47     # coming from taking a suitable root of a negative generator; this special
     48     # element is dealt with later). The shortlist grows at each step, so we
     49     # declare it before starting to loop over n and we build it incrementally.
     50     # The "special element" is of the form \zeta_{2^n}\sqrt{b}, which we encode
     51     # as (n,b). We use the value (1,1) (special element -1) to say that there
     52     # isno special element at this level.
     53     shortlist = []
     54     special_element = (1,1)
     55 
     56     # The integers M, giving the smallest cyclotomic field in which lies the
     57     # whole intersection with Q_\infty, also grows with n. As with the
     58     # shortlist, we declare it here and increase it appropriately at each step.
     59     M = 1
     60 
     61     for n in range( 1, N+1 ): # 1 \leq n \leq N
     62 
     63         # We add the new elements to the shortlist, modifying M if needed.
     64         # This is not done in case we are in the extra "fake" level (this case
     65         # dealt with immediately below).
     66         if n-1 < len(B):
     67             for g in B[n-1]:
     68                 # Case of negative g
     69                 if g < 0 and n > 1:
     70                     # Special element of the form \zeta_{2^{n+1}}\sqrt(b).
     71                     # It is contained in Q_{lcm(2^{n+1},cyc_emb(b))}, except in 
     72                     # case n=2 and b=2, in which case it is contained in Q_4.
     73                     # We store it as as pair ( n+1, b ).
     74                     special_element = ( n+1, abs(g)^(1/(2^(n-1))) )
     75                     M = lcm( M, special_embed( special_element ) )
     76                 else:
     77                     b = g^(1/(2^(n-1))) # b is 2-indivisible
     78                     shortlist.append( b )
     79                     M = lcm( M, cyc_embed(b) )
     80 
     81         # We add a root of an even power of the negative generator, as soon as
     82         # we are beyond its level.
     83         if d != -1 and n == d+2:        
     84             b = abs(B[d][0])^(1/2^d)
     85             shortlist.append( b ) 
     86             M = lcm( M, cyc_embed(b) )
     87 
     88         M = lcm(M,2^n)
     89             
     90         # Here we account for the extra 2^{n+1} root of unity (in case the
     91         # negative element has enough divisibility).
     92         if n <= d:
     93             M = lcm( M, 2^(n+1) )
     94 
     95         # We have to add -1 to the shortlist (see example [12,36]) and
     96         # remove it later.
     97         if n == 1 and d >= 1:
     98             shortlist.append(-1)
     99         if n > 1 and -1 in shortlist:
    100             shortlist.remove(-1)
    101 
    102         # For each divisor dM of M, compute the degree of the intersection of
    103         # Q(\sqrt{2^n}{G}) with Q_dM over Q_{2^n}. We need to compute the
    104         # number r of elements in the subgroup of G generated by the shortlist
    105         # that lie in Q_dM. This is going to be a power of 2. The sought degree
    106         # will be r, up to considering some special cases (described below).
    107         # 
    108         # This algorithm could be inefficient for groups of big rank. It can be
    109         # improved by precomputing subgroups of G/G^2 and the cyclotomic fields
    110         # containing them.
    111 
    112         aux = [] # Next line of ad_fail table
    113 
    114         for dM in divisors( M ):
    115             # We only care of the intersection with Q_dM if it contains the 2^n
    116             # roots of unity.
    117             if dM % (2^n) != 0:
    118                 continue
    119 
    120             # The following line gives, without repetitions, the list of m's
    121             # such that that some element in the subgroup of G generated by
    122             # the shortlist embeds in Q_m.
    123             S = [ product(s) for s in subsets( shortlist ) ]
    124             H = [ cyc_embed( s ) for s in S ] 
    125             r = len( [ b for b in H if dM % b == 0 ] )
    126             
    127             # We double n in case a new roots of unity enters the cyclotomic
    128             # due to the negative negerator. In case n=1, this is accounted
    129             # by having -1 in the shortlist.
    130             if n <= d and dM % (2^(n+1)) == 0 and n > 1:
    131                 r *= 2
    132 
    133             # We loose a factor of 2 if we have sqrt(2) in Q_8 or \zeta_8 2 in
    134             # Q_4.
    135             if 8 in H and dM % 8 == 0 and (n >= 3 or (n == 2 and n <= d)):
    136                 r = r/2
    137 
    138             # If we have a special element in this level, we consider it.
    139             # We have to consider all possible special elements arising from
    140             # multiplying the given element with the other elements of the
    141             # shortlist.
    142             # If any of them is of the form \zeta_8 2q for q a square, there is
    143             # nothing to do: in fact \zeta_8 2 embeds in Q_4, which coincides
    144             # with Q_{2^n} (n must be 2); so we would double the degree because
    145             # of the existence of the special element, butthen we would loose
    146             # another factor of 2 because of this.
    147             if special_element != (1,1) and special_element[0] == n+1:
    148                 nothing_to_do = False
    149                 intersecting_QdM = False
    150                 for s in S:
    151                     new_special = ( n+1, special_element[1] * s )
    152                     m = special_embed( new_special )
    153                     if n == 2 and m == 4: # \zeta_8 times 2 times square
    154                         nothing_to_do = True
    155                     if dM % m == 0:
    156                         intersecting_QdM = True
    157                 if intersecting_QdM and not nothing_to_do:
    158                     r *= 2
    159 
    160             aux.append( (dM,r) )
    161         
    162         ad_fail.append(aux)
    163 
    164     return ad_fail
    165 
    166 # Returns the smallest m such that \sqrt(b) is in the m-th cyclotomic field.
    167 def cyc_embed( b ):
    168     m = squarefree_part(b)
    169     if m%4 != 1:
    170         m *= 4
    171     return abs(m)
    172 
    173 # Computes the minimal cyclotomic field containing \zeta_{2^n}\sqrt(b),
    174 # where (n,b)=p
    175 def special_embed( p ):
    176     n = p[0]
    177     b = p[1]
    178     m = squarefree_part(b)
    179     if n == 3 and m % 2 == 0:
    180         return 4 * cyc_embed(m/2)
    181     else:
    182         return lcm( 2^n, cyc_embed(b) )
    183 
    184 # Computes the "l-adic failure", i.e. the ratio between l^{nr} and the degree
    185 # of Q_{l^m}(\sqrt[2^n](G)) over Q_{l^m} for all possible l.
    186 # Input: any basis for G
    187 # Returns a pair (L,T) where:
    188 #   - L is a list of primes l that do not have maximal l-adic part
    189 #   - T is an array of tables, where T[i][m][n] is the l-valuation of the
    190 #     degree of Q_{l^m,l^n} over Q_{l^m}, where l = L[i]
    191 def total_l_adic_failure( B ):
    192 
    193     bp = bad_primes( B )
    194     ret = ( bp, [] )
    195 
    196     for l in bp:
    197         ret[1].append( l_adic_failure( B, l ) )
    198 
    199     return ret
    200 
    201 # Computes the "bad primes", i.e. the ones for which the l-adic part is not
    202 # maximal. Used by l_adic_degree and total_l_adic_failure.
    203 def bad_primes( B ):
    204     M = exponent_matrix( B )
    205     (a,b) = M.dimensions()
    206     if a > b or M.rank() < a:
    207         print( "This is not a basis" )
    208         return
    209 
    210     # Compute which primes l divide all minors of the exponent matrix
    211     ms = M.minors( a )
    212     d = ms[0]
    213     for m in ms:
    214         d = gcd( d, m )
    215     bad_primes = list( prime_factors( d ) )
    216     if 2 not in bad_primes:
    217         bad_primes += [2] # 2 is always bad
    218     bad_primes.sort() # Ensures 2 is always first
    219 
    220     return bad_primes
    221    
    222 # Computes the l-adic failure for a specific l. Returns a "table" as described
    223 # above "total_l_adic_failure".
    224 # B is any basis for G.
    225 def l_adic_failure( B, l ):
    226 
    227     r = len(B)
    228     GB = make_good_basis( B, l )
    229     
    230     # Computes the parameters over Q4. For l odd, they are the same as over Q.
    231     p = parameters_Q4( GB, l )
    232     maxM = max( [ sum(x) for x in p ] )
    233     maxN = max( maxM, len(GB)-1 )
    234 
    235     maxM = max( 1, maxM )
    236     if l == 2:
    237         maxM = max( 2, maxM ) # For computing parameters over Q4
    238     maxN = max( 1, maxN )
    239 
    240     tabel = []
    241     max_failure = 0
    242     for m in range( 1, maxM+1 ):
    243         row = []
    244         for n in range( 1, max( m+1, maxN ) ):
    245             vl = -1
    246             if m >= n:
    247                 if l==2 and n==1 and m==1:
    248                     vl = compute_vl( p, 1, 2, r )
    249                     if adjust_sign( GB )[1] >= 1:
    250                         vl += 1
    251                 else:
    252                     vl = compute_vl( p, n, m, r )
    253                 failure = r*n - vl
    254             row.append(vl)
    255             max_failure = max( max_failure, failure )
    256         tabel.append((row,max_failure))
    257     return tabel
    258 
    259 # Returns a power of l that is the l-adic failure at M, N.
    260 # tablel must be the output table of l_adic_failure.
    261 def l_adic_failure_from_data( B, l, tablel, M, N ):
    262 
    263     m = valuation( M, l )
    264     n = valuation( N, l )
    265     if m < n:
    266         # inconsistent choice of M and N
    267         return 0
    268     if n == 0:
    269         return 1
    270     r = len(B)
    271 
    272     # Basically, the "dual" of what we do in l_adic_degree.
    273     if m > len(tablel):
    274         if n > len(tablel):
    275             return l^tablel[-1][1]
    276         else:
    277             return l^(r*n-tablel[-1][0][n-1])
    278     
    279     return l^(r*n-tablel[m-1][0][n-1])
    280 
    281 # Computes the l-divisibility parameters of G over Q4, given a good basis b
    282 # over Q for G. Returns a list of pairs (di,hi).
    283 # If l is odd it just uses the good basis given to compute the parameters.
    284 def parameters_Q4( gb, l ):
    285     # Converts from "good basis format" to simple list
    286     b = []
    287     for x in gb:
    288         b += x
    289     ret = []
    290     
    291     if l != 2:
    292         for i in range( len( gb ) ):
    293             for j in gb[i]:
    294                 ret.append( (i,0) )
    295         return ret
    296     else:
    297         R.<y> = PolynomialRing( QQ )
    298         pol = R(y^2+1)
    299         Q4.<eye> = NumberField( pol ) # I already use i for other things
    300 
    301         # Factorize basis elements over Q4 and so on.
    302         d = []
    303         B = []
    304         h = []
    305         ideals_list = set()
    306         M = [] # Exponent matrix of the Bi's
    307 
    308         # Pre-process to find all ideals appearing in the factorization and fix
    309         # a chosen generator for each of them. This is important in order to
    310         # compute the "sign" (h-parameter) of an element with respect to it Bi.
    311         for g in b:
    312             factorization_list = list( Q4.ideal(g).factor() )
    313             ideals_list |= set( [ x[0] for x in factorization_list ] )
    314         ideals_list = list( ideals_list )
    315         # Chooses a generator of each principal ideal in the list
    316         irreducibles_list = [ J.gens_reduced()[0] for J in ideals_list ]
    317 
    318         # Compute the Q4-parameters of the given basis b. Also computes the
    319         # exponent matrix of the Bi's 
    320         for g in b:
    321             factorization_list = list( Q4.ideal(g).factor() )
    322             exps = [ x[1] for x in factorization_list ]
    323             d.append( divisibility( exps, l ) )
    324             Bg = 1
    325             for j in range(len(factorization_list)):
    326                 a = 0
    327                 for i in range( len( ideals_list ) ):
    328                     if ideals_list[i] == factorization_list[j][0]:
    329                         a = irreducibles_list[i]
    330                         break
    331                 Bg *= a ^ (exps[j]/(l^d[-1]))
    332             B.append(Bg)
    333             u = g / (Bg^(l^d[-1]))
    334             if not u.is_unit():
    335                 print( "Error: g is not the right power of the computed Bg." )
    336                 print( "g:", g, ", Bg:", Bg, ", exponent:", l^d[-1] )
    337             if u == 1:
    338                 h.append( 0 )
    339             elif u == -1:
    340                 h.append( 1 )
    341             else:
    342                 h.append( 2 ) 
    343 
    344         # Make the exponent matrix M (for now as a list of rows)
    345         for g in B:
    346             row = [0] * len(ideals_list)
    347             for i in range(len(ideals_list)):
    348                 I = ideals_list[i]
    349                 ee = 1
    350                 while (I^ee).divides(g):
    351                     ee += 1
    352                 row[i] = ee-1
    353             M.append(row)
    354 
    355         # If the Bi's are not strongly independent, apply the algorithm (only
    356         # once) to produce a new basis. The new basis has maximal parameters.
    357         coeffs = find_combination( matrix(M), l )
    358         if coeffs != []:
    359 
    360             maxi = -1
    361             maxd = -1
    362             for i in range(len(d)):
    363                 if d[i] > maxd and coeffs[i] != 0:
    364                     maxd = d[i]
    365                     maxi = i
    366 
    367             x = [(a/coeffs[maxi]).lift() for a in coeffs] # Now a vector of int
    368 
    369             new_element = 1
    370             for i in range(len(d)):
    371                 new_element *= b[i]^( x[i] * l^(d[maxi]-d[i]) )
    372         
    373             b[maxi] = new_element
    374 
    375             # Compute new B, d and so on.
    376             factorization_list = list( Q4.ideal(b[maxi]).factor() )
    377             exps = [ x[1] for x in factorization_list ]
    378             d[maxi] = divisibility( exps, l )
    379             Bg = 1
    380 
    381             for j in range(len(factorization_list)):
    382                 a = 0
    383                 for i in range( len( ideals_list ) ):
    384                     if ideals_list[i] == factorization_list[j][0]:
    385                         a = irreducibles_list[i]
    386                         break
    387                 Bg *= a ^ (exps[j]/(l^d[maxi]))
    388             B[maxi] = Bg
    389             M[maxi] = [ x[1] for x in list( Q4.ideal(Bg).factor() ) ]
    390             u = b[maxi] / (Bg^(l^d[maxi]))
    391             if not u.is_unit():
    392                 print( "Error: new element is not the right power of B." )
    393                 print( "New el.:",b[maxi],", B:",Bg,", exponent:",l^d[maxi] )
    394             if u == 1:
    395                 h[maxi] = 0
    396             elif u == -1:
    397                 h[maxi] = 1
    398             else:
    399                 h[maxi] = 2 
    400             
    401         return [(d[i],h[i]) for i in range(len(d))]
    402 
    403 # Uses Theorem 18 to compute the degree of Kummer extensions.
    404 def compute_vl( p, n, m, r ):
    405     h = [ x[1] for x in p ]
    406     ni = [ min( n, x[0] ) for x in p ]
    407     M = max( m, max( [ h[i] + ni[i] for i in range( len( p ) ) ] ) )
    408     
    409     return M - m + r*n - sum( ni )
    410 
    411 # Given any basis b of a group G computes an l-good basis for G. This is done
    412 # using the algorithm outlined in the proof of Theorem 14 of (Debry-Perucca).
    413 def make_good_basis( b, l ):
    414     M = exponent_matrix( b )
    415     d = []
    416     B = []
    417     for i in range(len(b)):
    418         di = divisibility( M[i], l )
    419         d.append( di )
    420         B.append( abs(b[i])^(1/(l^di)) )
    421 
    422     # Computes the coeffiecients of a linear combination of the rows of M
    423     # that is zero modulo l. These coefficients are elements of F_l.
    424     coeffs = find_combination( exponent_matrix( B ), l )
    425 
    426     while coeffs != []:
    427 
    428         # Computes which basis element (with non-zero coefficient in the linear
    429         # combination above) has maximal divisibility.
    430         maxi = -1
    431         maxd = -1
    432         for i in range(len(d)):
    433             if d[i] > maxd and coeffs[i] != 0:
    434                 maxd = d[i]
    435                 maxi = i
    436 
    437         x = [ (a/coeffs[maxi]).lift() for a in coeffs ] # Now a vector of ints
    438 
    439         new_element = 1
    440         for i in range(len(d)):
    441             new_element *= b[i]^( x[i] * l^(d[maxi]-d[i]) )
    442         
    443         b[maxi] = new_element
    444         M = exponent_matrix( b )
    445         d[maxi] = divisibility( M[maxi], l )
    446         B[maxi] = abs(b[maxi])^(1/(l^d[maxi]))
    447         
    448         coeffs = find_combination( exponent_matrix( B ), l )
    449 
    450     GB = [[]]
    451     for i in range(len(d)):
    452         while( len(GB) <= d[i] ):
    453             GB.append([])
    454         GB[d[i]].append(b[i])
    455     return GB
    456 
    457 # Takes a good basis B and adjusts the sign of the elements so that there is at
    458 # most one negative generator (of positive divisibility). The input is a good
    459 # basis in the format returned by make_good_basis.
    460 # Returns a pair (B,d), where B is the updated basis and d is the divisibility
    461 # parameter of the only negative element remained.
    462 # The sign of the d=0 elements is just ignored in the other steps of the
    463 # algorithm, so we keep them negative.
    464 def adjust_sign( B ):
    465     neg = 0
    466     dret = -1
    467     for d in range( len(B)-1, 0, -1 ):
    468         for i in range(len(B[d])):
    469             if B[d][i] < 0:
    470                 if neg == 0:
    471                     neg = B[d][i]
    472                     dret = d
    473                 else:
    474                     B[d][i] *= neg
    475     return (B,dret)
    476 
    477 # Given the exponent matrix M of a list of rational numbers B, returns the
    478 # coefficients of a linear combination that is weakly l-divisible, or [] if
    479 # the B[i] are strongly l-independent.
    480 def find_combination( M, l ):
    481     M = M.change_ring( GF( l ) )
    482     if M.rank() != min( M.dimensions() ):
    483         return M.kernel().basis()[0]
    484     else:
    485         return []
    486 
    487 # A rational number must be given as a list of exponents of primes in its
    488 # factorization.
    489 # Returns the minimal l-valuation of the exponents.
    490 def divisibility( A, l ):
    491     if len(A) == 0:
    492         print( "Warning: computing the divisibility of a torsion element.", end="" )
    493         print( "Returning +Infinity." )
    494         return +Infinity
    495     return min( [ valuation( x, l ) for x in A ] )
    496 
    497 # For a given set of non-zero rationals B, computes the "exponent matrix"
    498 def exponent_matrix( B ):
    499     prime_list = set()
    500     for g in B:
    501         prime_list |= set( prime_factors( g ) )
    502     prime_list = list( prime_list )
    503     np = len( prime_list )
    504     rows = []
    505     for g in B:
    506         rowg = [0] * np
    507         for f in list( factor( g ) ):
    508             for i in range( np ):
    509                 if f[0] == prime_list[i]:
    510                     rowg[i] = f[1]
    511         rows.append( rowg )
    512     return matrix( rows )
    513 
    514 # Computing the exponent matrix, but keeping an extra column for the signs.
    515 # Returns a pair (M,L) where M is the modified exponent matrix and L is the
    516 # list of primes appearing in the factorization.
    517 def exponent_matrix_with_sign_and_primes( B ): 
    518     prime_list = set()
    519     for g in B:
    520         prime_list |= set( prime_factors( g ) )
    521     prime_list = list( prime_list )
    522     np = len( prime_list )
    523     rows = []
    524     for g in B:
    525         rowg = [0] * np
    526         for f in list( factor( g ) ):
    527             for i in range( np ):
    528                 if f[0] == prime_list[i]:
    529                     rowg[i] = f[1]
    530         s = 0
    531         if sgn(g) == -1:
    532             s = 1
    533         rowg.append( s )
    534         rows.append( rowg )
    535     return ( matrix( rows ), prime_list )
    536 
    537 # This is a wrapper function for total_kummer_failure( G, True ), see below.
    538 def TotalKummerFailure( G ):
    539     total_kummer_failure( G, True )
    540 
    541 # Input: any set of generators for a subgroup G of Q*.
    542 # If output=False, returns a 4-uple (t,MM,NN,D):
    543 # - t is a pair, where t[0] is the rank of G and t[1] is either True (if G has
    544 #   torsion) or False (is it does not).
    545 # - MM is the pair (M0,divisors(M0))
    546 # - NN is the pair (N0,divisors(N0))
    547 # - D is a table F, where F[j][i] is the ration between phi(m)n^r and the
    548 #   degree of Q_{m,n} (i.e., the Kummer failure at m,n, where m is the i-th
    549 #   divisor of M0 and n is the j-th divisor of N0 (in the computed list))
    550 #   computed for a torsion-free part of G.
    551 # If output is True, outputs this data in a human-readable way and does not
    552 # return any value.
    553 def total_kummer_failure( G, output ):
    554 
    555     # Computing a basis.
    556     (BM,BM_primes) = exponent_matrix_with_sign_and_primes( G )
    557     BM = BM.echelon_form()
    558     BM_primes.append(-1)
    559     B = []
    560     torsion = False
    561 
    562     for r in BM.rows():
    563         gr = product( [ BM_primes[i]^r[i] for i in range(len(r)) ] )
    564         if gr == -1:
    565             torsion = True
    566             break
    567         elif gr == 1:
    568             break
    569         else:
    570             B.append(gr)
    571 
    572     if len(B) == 0:
    573         print( "G is torsion. The extension is cyclotomic. Stopping." )
    574         return False
    575 
    576     r = len(B) # Rank of G
    577 
    578     # Compute l-adic data (straightforward)
    579     ( bad_primes, l_adic_failure_table ) = total_l_adic_failure( B )
    580 
    581     # Compute adelic data.
    582     (GB,d) = adjust_sign( make_good_basis( B, 2 ) )
    583     adelic_failure_table = adelic_failure_gb( GB, d )
    584     
    585     # Computing the bounds M0 and N0
    586     N0 = 1
    587     for i in range(len( bad_primes )):
    588         N0 *= bad_primes[i] ^ len( l_adic_failure_table[i][-1][0] )
    589     # Extra factors of 2 may come from the adelic failure
    590     N0 = lcm( N0, 2^len( adelic_failure_table ) )
    591     divs_N0 = divisors(N0)
    592 
    593     # greatest M appearing in the adelic failure table
    594     M0 = adelic_failure_table[-1][-1][0]
    595     divs_M0 = divisors(M0)
    596 
    597     # Failure Table
    598     FT = [ [ 1 for d1 in divisors(M0) ] for d2 in divisors(N0) ]
    599 
    600     # Takes the adelic failure from the the relative table
    601     for i in range(1,len(adelic_failure_table)+1):
    602         for pp in adelic_failure_table[i-1]:
    603             # Runs through all divisors of M0 and N0 that may have this
    604             # failure (the table is not "complete").
    605             for l in range(len(divs_M0)):
    606                 for j in range(len(divs_N0)):
    607                     dM = divs_M0[l]
    608                     dN = divs_N0[j]
    609                     if dN%(2^i)==0:
    610                         if dN%(2^(i+1))!=0 or i==len(adelic_failure_table):
    611                             if dM % pp[0] == 0:  
    612                                 FT[j][l] = lcm( FT[j][l], pp[1] )
    613 
    614     # Adding l-adic failure to the table
    615     for i in range( len( bad_primes ) ):
    616         l = bad_primes[i]
    617         for j in range(len(divs_N0)):
    618             dN = divs_N0[j]
    619             fl = l_adic_failure_from_data(B,l,l_adic_failure_table[i],dN,dN)
    620             for h in range(len(divs_M0)):
    621                 FT[j][h] *= fl
    622 
    623     ret = ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT )
    624 
    625     if output:
    626         print_total_table( ret )
    627         # Uncomment following line for case list description.
    628         #print_case_list( ret )
    629         return
    630     else:
    631         return ret
    632 
    633 # For the torsion tables, recall that when -1 is in G, the failure is defined
    634 # as the ration between 2^eN^r and the degree of Q_{M,N} over Q_M, where e=1
    635 # if N is even and e=0 otherwise.
    636 
    637 # Makes the failure table for the torsion case when M/N is even. In this case,
    638 # an entry of the table is doubled if the corresponding value of N (actually,
    639 # of gcd(N,N0) ) is even, and is kept the same otherwise.
    640 # The expected degree (over Q) 2^e * phi(M) * N^r, where e=1 if N is even and
    641 # e=0 otherwise.
    642 def torsion_table_even( data ):
    643 
    644     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
    645 
    646     new_FT = [ [ 1 for d1 in divs_M0 ] for d2 in divs_N0 ]
    647     for i in range(len(divs_N0)):
    648         for j in range(len(divs_M0)):
    649             if divs_N0[i] % 2 == 0:
    650                 new_FT[i][j] = 2 * FT[i][j]
    651 
    652     return new_FT
    653 
    654 # Makes the failure table for the torsion case when M/N is odd. In this case
    655 # the entry at (M,N) is taken from the torsion-free entry at (2M,N).
    656 # In other words, the expected degree (over Q) 2^e * phi(M) * N^r, where e=1 if
    657 # N is even and e=0 otherwise.
    658 def torsion_table_odd( data ):
    659 
    660     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
    661     
    662     new_FT = [ [ 1 for d1 in divs_M0 ] for d2 in divs_N0 ]
    663     new_data = ( ( r, False ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT )
    664     for i in range(len(divs_N0)):
    665         for j in range(len(divs_M0)):
    666             dN = divs_N0[i]
    667             dM = divs_M0[j]
    668             # The following lines just copy the failure from (2M,N) in the 
    669             # torsion-free case. It is easier to write it like this, although
    670             # it is not necessary in theory to compute the degree.
    671             exp_deg = euler_phi(2*dM) * dN^r
    672             deg = kummer_degree_from_total_table( 2*dM, dN, new_data ) 
    673             new_FT[i][j] = exp_deg / deg
    674     
    675     return new_FT
    676 
    677 def print_total_table( data ):
    678 
    679     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
    680     
    681     print( "M_0 =", M0 )
    682     print( "N_0 =", N0 )
    683     print( "" )
    684     print( "The following table shows the total failure of Kummer degrees", end="" )
    685     if torsion:
    686         print( "in\n case the quotient M/N is EVEN." )
    687     else:
    688         print( "." )
    689     print( "Columns correspond to values of M, rows to values of N" )
    690     print( "" )
    691     print( "The degree of the Kummer extension (M,N) can be extracted by taking" )
    692     print( "the value f (failure) of the entry at (gcd(N,N0),gcd(M,M0)) and" )
    693     print( "simply computing ed(M,N) / f, where ed(M,N) is the expected degree" )
    694     print( "of the Kummer extension." )
    695     if torsion:
    696         print( "In this case (-1 is in G), we have ed(M,N) = 2^e*phi(M)*N^r," )
    697         print( "where e=1 if N is even and e=0 if N is odd." )
    698         FT1 = torsion_table_even( data )
    699     else:
    700         print( "In this case (G is torsion-free) we have ed(M,N) = phi(M)*N^r," )
    701         FT1 = FT
    702     print( "where r is the rank of G." )
    703     print( "" )
    704 
    705     tt = [ ["","|"] + divs_M0 ]
    706     tt.append( "-" * (len(divs_M0)+2) )
    707     for i in range(len(divs_N0)):
    708         tt.append( [ divs_N0[i] ] + ["|"]  + FT1[i] )
    709     print( table(tt) )
    710     print( "" )
    711 
    712     if torsion:
    713         print( "The following table shows the total failure of Kummer degrees in" )
    714         print( "case the quotient M/N is ODD." )
    715         print( "This table can be read exactly as the first one." )
    716         print( "" )
    717 
    718         # A good strategy is the following:
    719         #   A little translation exercise: the failure at (M,N) in the torsion
    720         #   case is either the same as that for (2M,N) in the torsion-free case
    721         #   (if M is even) or its double (if M is odd).
    722         # However, due to problems in reading the table for bigger M, it is
    723         # easier to just compute the degree every time, and then deduce the
    724         # failure. This is not too inefficient, since we can use the data that
    725         # we have already computed via kummer_degree_from_total_table.
    726         new_FT = torsion_table_odd( data )
    727         # Printing the new table
    728         tt = [ ["","|"] + divs_M0 ]
    729         tt.append( "-" * (len(divs_M0)+2) )
    730         for i in range(len(divs_N0)):
    731             tt.append( [ divs_N0[i] ] + ["|"]  + new_FT[i] )
    732         print( table(tt) )
    733         print( "" )
    734 
    735 def print_case_list( data ):
    736 
    737     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
    738 
    739     FT_odd = torsion_table_odd( data )
    740     # FT1 is either FT or FT_even in the torsion case
    741     FT1 = FT
    742     if torsion:
    743         FT1 = torsion_table_even( data )
    744     pf = []
    745     for row in FT1:
    746         pf += row
    747     if torsion:
    748         for row in FT_odd:
    749             pf += row
    750     pf = list(set(pf))
    751     pf.sort()
    752     for f in pf:
    753         print( "Failure is", f, "if", end="" )
    754         if torsion:
    755             print( "M/N is EVEN and", end="" )
    756         print( "(gcd(M,M0),gcd(N,N0)) is one of the following:" )
    757         lijst = []
    758         for i in range(len(divs_N0)):
    759             for j in range(len(divs_M0)):
    760                 if FT1[i][j] == f:
    761                     lijst.append( ( divs_M0[j], divs_N0[i] ) )
    762         print( lijst )
    763         if torsion:
    764             print( "or if M/N is ODD and (gcd(M,M0)),gcd(N,N0)) is one of the", end="" )
    765             print( "following:" )
    766 
    767             lijst_odd = []
    768             for i in range(len(divs_N0)):
    769                 for j in range(len(divs_M0)):
    770                     if FT_odd[i][j] == f:
    771                         lijst_odd.append( ( divs_M0[j], divs_N0[i] ) )
    772             print( lijst_odd )
    773 
    774         print( "" )
    775 
    776 # Extracts a specific value of failure from the total table.
    777 def kummer_failure_from_total_table( M, N, data ):
    778     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
    779     FT1 = FT
    780     if torsion:
    781         if (M/N) % 2 == 0:
    782             FT1 = torsion_table_even( data )
    783         else:
    784             FT1 = torsion_table_odd( data )
    785 
    786     i = divs_N0.index( gcd( N, N0 ) )
    787     j = divs_M0.index( gcd( M, M0 ) )
    788     return FT1[i][j]
    789 
    790 # Computes the degree of the Kummer extension (M,N), by taking as input the
    791 # table computed by TotalKummerFailure.
    792 def kummer_degree_from_total_table( M, N, data ):
    793     
    794     ( ( r, torsion ), ( M0, divs_M0 ), ( N0, divs_N0 ), FT ) = data
    795 
    796     exp_deg = euler_phi(M) * N^r
    797     if torsion and N % 2 == 0:
    798         exp_deg *= 2
    799     return exp_deg / kummer_failure_from_total_table( M, N, data )
    800 
    801 # Given a set of generators for a finitely generated subgroup G of the
    802 # multiplicative group of Q, returns the degree of Q_{M,N} over Q.
    803 # M must be a multiple of N.
    804 def KummerDegree( G, M, N ):
    805     if M % N != 0:
    806         print( "M is not a multiple of N" )
    807         return -1
    808 
    809     data = total_kummer_failure(G,False)
    810     ((r,torsion),(M0,divs_M0),(N0,divs_N0),FT) = data
    811 
    812     exp_deg = euler_phi(M) * N^r
    813     
    814     if torsion and N % 2 == 0:
    815         exp_deg *= 2
    816     
    817     j = divs_M0.index(gcd(M,M0))
    818     i = divs_N0.index(gcd(N,N0))
    819 
    820     if torsion: 
    821         if (M/N)%2 == 0:
    822             failure = torsion_table_even( data )
    823         else:
    824             failure = torsion_table_odd( data )
    825     else:
    826         failure = FT
    827 
    828     return exp_deg / failure[i][j]