/************************************************************************************** Input: K real quadratic field l small prime number (l=3,5,7) E l-th cyclotomic field rcgp ray class group with conductor supported in P subgrp subgroup of rcgp corresponding to Kp Description: Verification of ETNC as described in the paper "Computation of Stark-Tamagawa units" The implementation only works under the following HYPOTHESIS: (i) h_K = 1, h_E = 1, (ii) [K:Q] = 2 (iii) the ramified primes generate the ideal class group of Kp Output: vector v of length 20 v[1] Kp as output by bnfinit v[2] S-units of Kp, where S = S_infty + S_ram v[3] S-class number (which must be 1 to obtain a correct result) v[4] vector of S-units generating the lattice Eps_S v[5] index of Eps_S in the S-units mod torsion v[6] vector of determinants of the regulator matrices v[7] vector of L-values v[8] vector of the complex numbers a_chi v[9] 2-component vector v[9][1] embedding of K into Kp v[9][2] primes of Kp above P v[10] representative of rcgp mod subgroup v[11] isomorphism iso: Kp --> Kp such that the Artin symbol of v[10] equals iso v[12] representation (over Z) induced by v[2] v[13] 2-component vector of Galois orbit representatives of the infinite places of Kp v[14] [eps_1, eps_2], the units which generate e_1*U_S over e_1*ZG v[15] [eta_1, ..., eta_r], units in K such that has index in U_S prime to l and such that (eta_i, Kp_j/K_j) = delta_ij v[16] polynomial prod_{chi != chi_0} (X - a_chi), already rounded v[17] maximal rounding error v[18] 2-component vector [a_chi_0, a_chi] as an element in Q + E v[19] not used v[20] not used *****************************************************************************************/ etnc( K, l, E, rcgp, subgrp) = { local( h, Kp, v, embed, PKp, p, KpSunits, xi, i, j, Sunits, P, result, r ); result = vector(20, i, 0); moduleFact = idealfactor(K, rcgp[2][1]); P = vector(matsize(moduleFact)[1], k, moduleFact[k,1]); print( "Compute Kp (the cyclic extension of K of degree l)"); h = bnrstark( rcgp, subgrp, 2 ); Kp = bnfinit( h ); result[1] = Kp; print("Compute embedding of K into Kp and prime ideals above P"); p = vector(length(P), i, P[i][1]); v = findPKp(K, Kp, P, l); embed = v[1]; PKp = v[2]; result[9] = v; print("Compute S-units"); KpSunits = bnfsunit( Kp, PKp ); xi = vector(length(p), j, Mod(KpSunits[1][j], Kp.pol)); /* Dies sind die zusaetzlichen S-Einheiten */ Sunits = vector(2*l-1+length(p), i, 0.0); for (i=1, 2*l-1, Sunits[i] = Mod(Kp.fu[i], Kp.pol)); for (j=1, length(p), Sunits[2*l-1+j] = xi[j]); result[2] = Sunits; print("Compute S-class number"); hS = SClassNumber( Kp, PKp ); if( gcd( hS, l ) != 1, print( " S-Klassengruppe hat Ordnung\n nicht teilerfremd zu l " ); return( 0 ); ); result[3] = hS; print( "Compute the Galois group " ); allisos = CompGal(Kp, K, embed); vertreter = findvertreter( rcgp, subgrp, l ); iso = findiso( K, Kp, rcgp, l, vertreter, allisos, embed ); DarSigma = RepOfG( Kp, iso, KpSunits ); result[10] = vertreter; result[11] = iso; result[12] = DarSigma; print("Compute archimedian places"); alpha = FindInfinitePlaces( Kp, l, iso ); result[13] = alpha; print( "Compute eps_1 and eps_2 (archimedian units)" ); archunits = ArchUnits(Kp, KpSunits, E, DarSigma, l); result[14] = archunits; print("Compute units eta_1, ..., eta_r (S-units in K)"); Eta = FindSUnitsOfK( K, Kp, embed, KpSunits, l, DarSigma, archunits, p, P ); Etaalg = unitbasistoalg( Eta, Sunits ); for(i=1, length(Etaalg), Etaalg[i]=setdown(K, Kp, Etaalg[i], embed )); if (length(P) != 1, C = FindC(K, l, Etaalg, moduleFact, rcgp, vertreter, subgrp): Etacorr = CorrUnitsOfK( Eta, C ); Etacorralg = unitbasistoalg( Etacorr, Sunits ) , Etacorr = Eta; Etacorralg = Etaalg; ); result[15] = vector(length(Etacorralg), i, setdown(K, Kp, Etacorralg[i], embed )); EpsMat = EpsBasisRepresentation( archunits, Etacorr, DarSigma, l ); Eps = concat( unitbasistoalg( archunits, Sunits ), Etacorralg ); result[4] = Eps; result[5] = matdet(EpsMat); print("Compute regulator matrices"); M = RegMat( Kp, Eps, alpha, PKp, l, DarSigma ); result[6] = vector(length(M), i, matdet(M[i])); print( "Compute L-values " ); Lval = bnrL1( rcgp, subgrp, 6 ); result[7] = vector(length(Lval), i, if (i==1, Lval[l][2][2], Lval[i- 1][2][2])); print("Compute a_chi"); z = exp( 2*Pi*I / l ); r = length( p ) - 1; B = vector( l, i, 0 ); for( i = 1, length( B ), if( i == 1, B[i] = matdet( M[i] )/( l^r * Lval[l][2][2] ), B[i] = ( z^( i - 1 )- 1)^r * matdet( M[i] ) / Lval[l - i + 1][2][2]; ); ); result[8] = B; print("Verification of ETNC"); RM = 1; for( i = 2, length( B ), RM = RM*( x - B[i] ); ); err1 = absdiff( RM ); /* Das groesste Rundungsfehler koeffizientenweise */ RRM = round( real( RM ) ); result[16] = RRM; FRM = factor( RRM )[ 1, 1 ]; iota = nfisincl( FRM, E.pol )[1]; b = Sigma( x, E.pol, iota ); PE = idealprimedec( E, l )[1]; IV = idealval( E, b - round( real( B[1] ) ), PE ); b0 = round( real(B[1])); err0 = abs( B[1] - b0 ); result[17] = max(err1, err0); result[18] = [b0, b]; if (max(err1, err0) < 10^(-(floor(precision(B[1])/2))) && abs(b0) == 2*abs(matdet(EpsMat)) && IV > 0, print("ETNC is valid"); return([1, result]); , print("COUNTER EXAMPLE (you should doubt the correctness of the computation)"); return([0, result]); ); } /****************************************************************************************/ /* *************** FUNCTIONS RELATED TO THE SPECIFIC PROBLEM ************************** */ /****************************************************************************************/ /* CAUTION: h_K = 1 is needed */ findPKp(K, Kp, P, l) = { local( p, i, i0, a, allembeds, Q, PKp, j, beta, embed ); \\ p = vector(length(P), i, idealnorm(K, P[i])); p = vector(length(P), i, P[i][1]); allembeds = nfisincl( K.pol, Kp.pol ); Q = vector(length(p), i, idealprimedec(Kp, p[i])); PKp = vector(length(p), j, -1); i0 = -1; for (i=1, length(p), for (j=1, length(Q[i]), if (Q[i][j][3] == 1, i0 = i; break); ); if (i0 > 0, break); ); if (i0 < 0, i0 = 1); a = nfbasistoalg(K, bnfisprincipal(K, P[i0])[2]); a = subst(lift(a), y, x); for (i=1, length(Q[i0]), if (Q[i0][i][3]/P[i0][3] == l, PKp[i0] = Q[i0][i]; break)); for (i=1, length(allembeds), beta = Sigma(a, Kp.pol, allembeds[i]); if (idealval(Kp, beta, PKp[i0]) != 0, embed = allembeds[i]; break); ); for (i=1, length(P), if ( i != i0, a = nfbasistoalg( K, bnfisprincipal(K, P[i])[2] ); a = subst(lift(a), y, x); beta = Sigma(a, Kp.pol, embed); for (j=1, length(Q[i]), if (idealval(Kp, beta, Q[i][j]) != 0, PKp[i] = Q[i][j]; break); ); ); ); return([embed, PKp]); } CompGal(L, K, embed ) = { local( alliso, G, i ); alliso = nfgaloisconj(L); G = []; for (i=1, length(alliso), if (Mod(embed, L.pol) == Sigma(embed, L.pol, alliso[i]), G = concat(G, [alliso[i]])); ); return(G); } /* This function computes an isomorphism iso: bnfext --> bnfext such that the Artin-Symbol of representative equals iso */ findiso( bnf, bnfext, rnf, ordG, representative, autbnfext, embed ) = { local( i, j, k, tmp, factorlist, puffer ); local( relnorm, q, Q, A, B, C, beta, Qlist ); factorlist = idealfactor( bnf, representative ); puffer = vector( matsize( factorlist )[1], i, 0 ); for( k = 1, matsize( factorlist )[1], tmp = factorlist[k, 1]; relnorm = idealnorm( bnf, tmp ); q = factor( relnorm )[1, 1]; beta = nfbasistoalg(bnf, bnfisprincipal(bnf, tmp)[2]); beta = subst(lift(beta), y, x); beta = Sigma(beta, bnfext.pol, embed); Qlist = idealprimedec( bnfext, q ); i=1; while(idealval(Kp, beta, Qlist[i]) == 0, i++); Q = Qlist[i]; i = 0; while( !( i == length(autbnfext)), i = i + 1; j = 0; while( !( j == length( bnfext.zk ) ), j = j + 1; A = Sigma( bnfext.zk[j], bnfext.pol, autbnfext[i] ); A = rest_mod( bnfext, A, Q ); B = Mod( bnfext.zk[j], bnfext.pol ); B = power( bnfext, B, relnorm, Q ); C = rest_mod( bnfext, A - B, Q ); if( !( C == 0 ), j = length( bnfext.zk ); ); ); if( C == 0, puffer[k] = Mod( autbnfext[i], bnfext.pol ); i = length(autbnfext); ); ); ); for( i = 1, length( puffer ), /* Berechne die Potenzen von Automorphismen */ upper = ( factorlist[i, 2]%ordG ); tmp = puffer[i]; puffer[i] = Mod( x, bnfext.pol ); for( j = 1, upper, puffer[i] = subst( lift( puffer[i] ), x, tmp ); ); ); tmp = puffer[1]; for( i = 2, length( puffer ), /* Berechne Produkt von Automorphismen */ tmp = subst( lift( tmp ), x, puffer[i] ); ); tmp = lift( tmp ); return(tmp); } /* Choose infinite places such that their orbits form the full set of infinite places */ FindInfinitePlaces( bnf, ordG, aut ) = { local( flag, tolerance, bahn, tmp, alpha, kandidat, i, j ); bahn = List( [] ); tolerance = 0.001; for ( i = 1, ordG, tmp = conjvec( Sigmapow( Mod( x, bnf.pol ), bnf.pol, aut, i ) ); bahn = concat( bahn, List( [ tmp[1] ] ) ); ); alpha = tmp[1]; for ( i = 1, length( tmp ), kandidat = tmp[i]; j = 1; flag = 0; while( j <= length( bahn ), if( abs( kandidat - bahn[j] ) < tolerance, flag = 1; ); j = j + 1; ); if( flag == 0, return( [alpha, kandidat] ); ); ); } /* Compute the representation of G = induced by the S-units */ RepOfG( bnf, iso, bnfSunits ) = { local( i, j, A, tmp, marke, size, fullSunits ); fullSunits = concat( bnfunit( bnf )[1], bnfSunits[1] ); size = length( fullSunits ); A = matrix( size, size, i, j, 0 ); for( i = 1, size, tmp = Sigma( fullSunits[i], bnf.pol, iso ); tmp = bnfissunit( bnf, bnfSunits, tmp ); marke = length( tmp ) - length( bnfSunits[1] ); for( j = 1, length( tmp ), if( j < marke, /* Torsionskomponente */ A[j, i] = tmp[j], /* wird weggelassen */ if( !( j == marke ), A[j - 1, i] = tmp[j]; ); ); ); ); return( A ); } ArchUnits(Kp, KpSunits, E, DarSigma, l) = { local( FU, RA, AU); FU = QuasiMinkowski( Kp, KpSunits, DarSigma, l ); RA = QuasiMinkowskiRepresentation( Kp, E, FU, DarSigma, l ); AU = CompArchEps( Kp, E, RA, FU, DarSigma, l ); return(AU); } /* Compute S-units xi_1, xi_2 such that they span over e_1Q[G] the lattice e_1*U_S. */ QuasiMinkowski( bnf, bnfSunits, D, ordG ) = { local( i, j, e, upper, tmp, idemnew ); local( marke, xi1, xi2, X1, X2, U ); local( idem0, id, idem1 ); idem0 = [ vector( l, i, 1 ), l]; id = [ vector( l, i, 0 ), 1]; id[1][1] = 1; idem1 = Add( id, [ ( -1 )*idem0[1], l ] ); upper = 2*ordG + length( bnfSunits[1] ) - 1; idemnew = idem1; idemnew[2] = 1; /* Bei der Berechnung von Spann verzichte auf den Nenner */ i = 1; while( i <= upper, tmp = [ vector( upper, e, 0 ), 1 ]; tmp[1][i] = 1; xi1 = tmp; xi1 = groupringaction( idem1, xi1, D ); tmp = groupringaction( idemnew, tmp, D ); if( testnull( tmp ), i = i + 1, marke = i; i = upper + 1; ); ); X1 = unitspan( xi1, D, ordG ); U = mathnf( X1, 1 )[2]; X1 = (1/xi1[2])*X1; /* Der Nenner wird beruecksichtigt */ i = marke + 1; while( i <= upper, tmp = [ vector( upper, e, 0 ), 1 ]; tmp[1][i] = 1; xi2 = tmp; /*???*/xi2 = groupringaction( idem1, xi2, D );/*???*/ tmp = groupringaction( idemnew, tmp, D ); tmp[1] = tmp[1]*U; j = 1; /* Test "xi2 im Spann von xi1"*/ while( j <= ( upper - ordG + 1 ) && tmp[1][j] == 0, j = j + 1; ); if( j <= upper - ordG + 1, i = upper + 1; ); i = i + 1; ); X2 = unitspan( xi2, D, ordG ); X2 = (1/xi2[2])*X2; /*Der Nenner wird beruecksichtigt*/ return([[xi1, X1], [xi2, X2]]); } /* Computes the matrix A with coefficients in E = Q(zeta_l) such that e_1*(basis of S-units) = (xi_1, xi_2)*A */ QuasiMinkowskiRepresentation( bnf, zyknf, T, D, ordG ) = { local( evec, tmp, upper, i, j, A, X, v ); local( idem0, id, idem1 ); idem0 = [ vector( l, i, 1 ), l]; id = [ vector( l, i, 0 ), 1]; id[1][1] = 1; idem1 = Add( id, [ ( -1 )*idem0[1], l ] ); upper = matsize( D )[1]; X = matrix( upper, 2*ordG - 2, i, j, if( j <= ordG - 1, T[1][2][j, i], T[2][2][j - ordG + 1, i] ) ); evec = vector( upper, i, vector( upper, j, 0 ) ); for( i = 1, upper, tmp = [vector( upper, j, 0 ), 1]; tmp[1][i] = 1; tmp = groupringaction( idem1, tmp, D ); evec[i] = (1/tmp[2])*tmp[1]; ); v = vectorv( upper, i, matinverseimage( X, evec[i]~ ) ); A = matrix( 2, upper, i, j, 0 ); for( j = 1, upper, A[1,j] = nfbasistoalg( zyknf, vectorv( ordG - 1, i, v[j][i] ) ); A[2,j] = nfbasistoalg( zyknf, vectorv( ordG - 1, i, v[j][i + ordG - 1] ) ); ); return( A ); } /* Construction of the S-units eps_1, eps_2 corresponding to the archemedian places. */ CompArchEps( bnf, zyknf, A, X, D, ordG ) = { local( i, j, z, s, w1, w2, U, H, HU, idem0 ); local( tmp1, tmp2, x1, x2, Q, vec, eps1, eps2 ); idem0 = [vector( ordG, i, 1 ), ordG]; z = matsize( A )[1]; s = matsize( A )[2]; A = ordG*A; HU = mathnf_allg( zyknf, A ); H = matrix( z, z, i, j, [concat( nfalgtobasis( zyknf, HU[1][i, s - 2 + j] )~, [0] ), ordG]); w1 = groupringaction( H[1,1], X[1][1], D ); tmp1 = groupringaction( H[1,2], X[1][1], D ); tmp2 = groupringaction( H[2,2], X[2][1], D ); w2 = add( tmp1, tmp2 ); Q = idealprimedec( zyknf, ordG )[1]; U = matrix( s, z, i, j, rest_mod( zyknf, HU[2][i, s - 2 + j], Q ) ); U = matrix( s, z, i, j, polcoeff( lift( U[i,j] ), 0 ) ); vec = vector( s, i, [vector( s, j, 0 ), 1] ); for( i = 1, length( vec ), vec[i][1][i] = 1; vec[i] = groupringaction( idem0, vec[i], D ); ); x1 = [vector( s, i, 0 ), 1]; x2 = [vector( s, i, 0 ), 1]; for( i = 1, length( vec ), tmp1 = scalmul( U[i,1], vec[i]); x1 = add( x1, tmp1 ); ); for( i = 1, length( vec ), tmp1 = scalmul( U[i,2], vec[i]); x2 = add( x2, tmp1 ); ); eps1 = add( w1, x1 ); eps2 = add( w2, x2 ); if( abs( eps1[2] ) != 1 || abs( eps2[2] ) != 1, print( " WARNING: die Einheiten haben Nenner! " ); ); return( [eps1, eps2] ); } /* Computation of the S-units eta_1, ..., eta_r of K such that norm(eps_1), norm(eps_2), eta_1, ..., eta_r generate the S-units of K up to index prime to l */ FindSUnitsOfK( mybnf, bnfext, embeds, bnfextSunits, ordG, D, eps, p, Pp ) = { local( i, j, marke, tetavec, tmpvec, tmp ); local( n, m, e0AU, b, Nu, a, A, H, ind ); local( idem0, id, idem1 ); idem0 = [ vector( l, i, 1 ), l]; id = [ vector( l, i, 0 ), 1]; id[1][1] = 1; idem1 = Add( id, [ ( -1 )*idem0[1], l ] ); marke = length( bnfextSunits[1] ); tmpvec = vector( marke + 1, i, 0 ); tmpvec[1] = Sigma( subst( mybnf.fu[1], y, x ), bnfext.pol, embeds ); for( i = 2, length( tmpvec ), tmp = nfbasistoalg( mybnf, bnfisprincipal( mybnf, Pp[i - 1] )[2] ); tmpvec[i] = Sigma( subst( lift( tmp ), y, x ), bnfext.pol, embeds ); ); tetavec = vector( length( tmpvec ), i, vector( 2*ordG + marke - 1, j, 0 ) ); for( i = 1, length( tetavec ), tmp = bnfissunit( bnfext, bnfextSunits, tmpvec[i] ); for( j = 1, length( tmp ), if( j < length( tmp ) - marke, /* Torsionskomponente wird weggelassen*/ tetavec[i][j] = tmp[j], if( j > length( tmp ) - marke, tetavec[i][j - 1] = tmp[j]; ); ); ); ); n = 2*ordG + length( bnfextSunits[1] ) - 1; m = length( tetavec ); e0AU = vector( length( eps ), i, groupringaction(idem0, eps[i], D)); b = vector( length( eps ), i, vectorv( n, j, e0AU[i][1][j] / e0AU[i][2] ) ); Nu = matrix( n, m, i, j, tetavec[j][i] ); a = vector( length( eps ), i, matinverseimage( Nu, b[i] ) ); A = matrix( m, 2, i, j, ordG*a[j][i] ); H = mathnf( A ); tmp = findkomp( H, ordG ); /* Wahl vom Minor mit der kleinsten Determinante */ tmpvec = vector( length( tmp ), i ); for( i = 1, length( tmp ), tmpvec[i] = [tetavec[tmp[i]], 1]; ); return( tmpvec ); } /* Compute the matrix C needed to adapt the S-units eta_1, ..., eta_r Explicitly, C = D^-1, where (teta_i, L_j/K_j) = g_0^{d_{ij}}. CAUTION: it is absolutely vital that h_K = 1 */ OldFindC(K, l, teta, P, rcgp, gen, subgrp)= { local( r, C, vgen, a, i, j, k, e, p, v, NullVec ); r = length(teta); C = matrix(r, r); vgen = bnrisprincipal(rcgp, gen)[1]; a = vector(r+1, i, 0); NullVec = vectorv(length(rcgp.clgp[2]), i, 0); for (i=1, r, for (j=1, r, e = idealval(K, teta[i], P[j+1]); p = nfbasistoalg(K, bnfisprincipal(K, P[j+1])[2]); for (k=0, r, if (k != j, a[k+1] = p^e, a[k+1] = InverseModId(K, teta[i] / (p^e), P[j+1])); ); v = china(K, a, P); v = bnrisprincipal(rcgp, v)[1]; for (c=0, l-1, if (VectorModHNF(v - c*vgen, subgrp) == NullVec, C[i, j] = Mod(c, l); break; ); ); ); ); return( lift(C^(-1)) ); } /* Compute the matrix C needed to adapt the S-units eta_1, ..., eta_r Explicitly, C = D^-1, where (teta_i, L_j/K_j) = g_0^{d_{ij}}. CAUTION: it is absolutely vital that h_K = 1 */ FindC(K, l, teta, moduleFact, rcgp, gen, subgrp)= { local( r, C, vgen, a, i, j, k, e, p, v, NullVec, IdList, P ); r = length(teta); C = matrix(r, r); vgen = bnrisprincipal(rcgp, gen)[1]; a = vector(r+1, i, 0); NullVec = vectorv(length(rcgp.clgp[2]), i, 0); P = vector(matsize(moduleFact)[1], k, moduleFact[k,1]); for (i=1, r, for (j=1, r, e = idealval(K, teta[i], P[j+1]); p = nfbasistoalg(K, bnfisprincipal(K, P[j+1])[2]); IdList = vector(matsize(moduleFact)[1], k, idealpow(K, moduleFact[k,1], moduleFact[k,2])); for (k=0, r, if (k != j, a[k+1] = p^e, a[k+1] = InverseModId(K, teta[i] / (p^e), IdList[j+1])); ); v = china(K, a, IdList); v = bnrisprincipal(rcgp, v)[1]; for (c=0, l-1, if (VectorModHNF(v - c*vgen, subgrp) == NullVec, C[i, j] = Mod(c, l); break; ); ); ); ); return( lift(C^(-1)) ); } /* Replace Eta by Eta * C */ CorrUnitsOfK( Eta, C ) = { local( i, j, k, len, size, tetatmp, puffer ); tetatmp = Eta; len = length( Eta[1][1] ); size = matsize( C )[1]; for( i = 1, size, puffer = vector( len, j, 0 ); for( j = 1, len, for( k = 1, size, puffer[j] = puffer[j] + C[i, k]*Eta[k][1][j]; ); ); tetatmp[i][1] = puffer; ); return( tetatmp ); } /* Compute the basis representation of Eps = _ZG with respect to the chosen basis of S-units. */ EpsBasisRepresentation( A, teta, D, ordG ) = { local( e, i, j, k, size, tmp, EPS ); size = matsize( D )[1]; EPS = matrix( size, size, i, j, 0 ); j = 0; for( e = 1, length( A ), /* Span von Einheiten zu archimedischen Stellen */ tmp = mattranspose( A[e][1] ); for( i = 1, ordG, tmp = D*tmp; j = j + 1; for( k = 1, size, EPS[k, j] = tmp[k]; ); ); ); i = size - length( teta ); for( e = 1, length( teta ), i = i + 1; for( j = 1, size, EPS[ j, i ] = teta[e][1][j]; /* Die Einheiten vom Grundkoerper */ ); ); return( EPS ); } /* Berechnung von logarithmischen Summen fuer verschiedene Charaktere von Galoisgruppe G. 0 <= flag < ordG. flag = 0 entspricht dem trivialen Charakter von G. Sunits ist ein Vek- tor der S-Einheiten in Kp, gegeben als polmods, embed ist die ausgewaehlte reelle Einbettung des Grundkoerpers. eps ist die in Kp liegende Einheit, deren logsum zu bestimmen ist. DarSigma ist die Darstellung von Galoisautomorphismus auf S-Einheiten in Kp. eps gegeben als polmod. bnfext entspricht Kp, S = Menge der nichtarchimedischen Primstellen. */ logsum( bnfext, S, DarSigma, ordG, embed, flag, eps ) = { local( i, j, k, var, len, tmp, prodtmp ); local( xi, wert, unitsval, epsnew, Sunits ); Sunits = concat( bnfext.fu, bnfsunit( bnfext, S )[1] ); xi = exp( 2*Pi*I/ordG ); var = variable( Sunits[1] ); len = length( Sunits ); unitsval = vector( len, i, 0 ); epsnew = vector( len, i, 0 ); tmp = bnfissunit( bnfext, bnfsunit( bnfext, S ), eps ); for( i = 1, len, unitsval[i] = subst( Sunits[i], var, embed ); ); for( i = 1, length( tmp ), if( i < 2*ordG, /* Nehme die Z-Koordinaten und */ epsnew[i] = tmp[i], /* lasse dabei die Torsion weg */ if( i != 2*ordG, epsnew[i - 1] = tmp[i]; ); ); ); wert = 0.0; if( flag != 0, k = ordG - flag, k = flag; ); tmp = mattranspose( epsnew ); /* Koordinatenkomponente */ /*print("Vor der eigentlichen Berechnung!\n");*/ for( i = 0, ordG - 1, if( i == 0, prodtmp = 1.0; for( j = 1, len, prodtmp = prodtmp*( unitsval[j]^tmp[j] ); ); wert = wert + log( abs( prodtmp ) )*xi^( k*i ), tmp = DarSigma*tmp; prodtmp = 1.0; for( j = 1, len, prodtmp = prodtmp*( unitsval[j]^tmp[j] ); ); wert = wert + log( abs( prodtmp ) )*xi^( k*i ); ); ); return( wert ); } /* Computes the regulator matrices R_\chi(\calE_S), where chi is determined by g_0 --> exp(2*Pi*I / l). */ OldRegMat( bnfext, eps, realembed, prim, ramprim, ordG, DarSigma ) = { local( M, M1, R, m, i, j, k, tmp ); M = List( [] ); M1 = matrix( m = length( ramprim ) + 1, m, i, j, 0.0 ); k = 0; for( i = 1, m, for( j = 1, 2, M1[i, j] = -logsum( bnfext, ramprim, DarSigma, ordG, realembed[j], k, eps[i] ); ); ); for( i = 1, m, for( j = 3, m, tmp = ( -idealval( bnfext, eps[i], ramprim[j - 1] ) ); M1[i, j] = -log( prim[j - 1]^tmp ); ); ); M = concat( M, List( [M1] ) ); R = matrix( 2, 2, i, j, 0.0 ); for( k = 1, ordG - 1, for( i = 1, 2, for( j = 1, 2, R[i, j] = -logsum( bnfext, ramprim, DarSigma, ordG, realembed[j], k, eps[i] ); ); ); M = concat( M, List( [R] ) ); ); return( M ); } RegMat( bnfext, eps, realembed, ramprim, ordG, DarSigma ) = { local( M, M1, R, m, i, j, k, tmp, Nv ); M = List( [] ); M1 = matrix( m = length( ramprim ) + 1, m, i, j, 0.0 ); k = 0; for( i = 1, m, for( j = 1, 2, M1[i, j] = -logsum( bnfext, ramprim, DarSigma, ordG, realembed[j], k, eps[i] ); ); ); for( i = 1, m, for( j = 3, m, tmp = ( -idealval( bnfext, eps[i], ramprim[j - 1] ) ); Nv = idealnorm(bnfext, ramprim[j - 1]); M1[i, j] = -log( Nv^tmp ); ); ); M = concat( M, List( [M1] ) ); R = matrix( 2, 2, i, j, 0.0 ); for( k = 1, ordG - 1, for( i = 1, 2, for( j = 1, 2, R[i, j] = -logsum( bnfext, ramprim, DarSigma, ordG, realembed[j], k, eps[i] ); ); ); M = concat( M, List( [R] ) ); ); return( M ); } /* Bestimmung der Untergruppen im Strahl mit dem Index ordG */ findsubgroup( rcgp, ordG ) = { local( i, tmp, sgl ); tmp = subgrouplist( rcgp, ordG, 0 ); sgl = List( [] ); for( i = 1, length( tmp ), if( matdet( tmp[i] ) == ordG, sgl = concat( sgl, List( [ tmp[i] ] ) ); ); ); return( sgl ); } /**************************************************************************/ /* *************** GROUPRING RELATED FUNCTIONS ************************** */ /**************************************************************************/ {Add(a, b, c, k, d, i, j) = c = [0, 0]; k = lcm(a[2], b[2]); if (type(a[1]) == "t_VEC", c[1] = k/a[2] * a[1] + k/b[2] * b[1]; c[2] = k; d = gcd(c[2], VecGcd(c[1])); for (i=1, length(a[1]), c[1][i] = c[1][i] / d); c[2] = c[2] / d; return(c); , c[1] = a[1]^(k/a[2]) * b[1]^(k/b[2]); c[2] = k; ); return(c); } /* Computation of the action of a group ring element alpha on the S-unit unit. unit is given by [eps, denom], where eps is a S-unit in basis representation and denom a denomnator in Z. alpha is a groupring element of the form [[a_0, ..., a_{l-1}], d], meaning \sum_{i=0}^{l-1}a_ig_0^i \in QG. */ groupringaction( alpha, unit, D ) = { local( i, lenal, lenun, tmp, size, B, T, ggT ); lenal = length( alpha[1] ); lenun = length( unit[1] ); T = D^0; B = alpha[1][1]*T; for( i = 2, lenal, T = T*D; B = B + alpha[1][i]*T; ); unit[1] = mattranspose( B*mattranspose( unit[1] ) ); unit[2] = unit[2]*alpha[2]; ggT = unit[2]; for( i = 1, lenun, ggT = gcd( unit[1][i], ggT ); ); for( i = 1, lenun, unit[1][i] = unit[1][i]/ggT; ); unit[2] = unit[2]/ggT; return( unit ); } /* ********************* MISCELLANOUS ************************************* */ SClassNumber( bnf, S ) = { local( A, i, j, n, m, tmp ); n = length( bnf.clgp.cyc ); m = n + length( S ); A = matrix( n, m, i, j, 0 ); for( i = 1, length( S ), tmp = bnfisprincipal( bnf, S[i] )[1]; for( j = 1, length( tmp ), A[j, i] = tmp[j]; ); ); j = 1; for( i = length( S ) + 1, m, A[j, i] = bnf.clgp.cyc[j]; j = j + 1; ); A = mathnf( A ); return( matdet( A ) ); /* Ordnung des Quotienten( cl(bnf)/cl(S) )*/ } /* Returns an ideal representing the ray class group mod subgroup if the quotient is cyclic of order ordG. */ findvertreter( rcgp, subgroup, ordG ) = { local( j ); j = 1; while( subgroup[j, j] != ordG, j = j + 1; ); return( rcgp.clgp[3][j] ); } /* Transorm the S-unit unit given in basis representation with respect to SU into an algebraic number. */ unitbasistoalg( unit, SU ) = { local( i, j, puffer, len, tmp ); len = length( unit ); puffer = vector( len, i, 0 ); for( i = 1, len, tmp = 1; for( j = 1, length( SU ), tmp = tmp*( SU[j]^unit[i][1][j] ); ); puffer[i] = tmp; ); return( puffer ); } /* If alpha \in bnf is given as an element of bnfext, this function computes alpha as an elemnt of bnf. CAUTION: only for [bnfext : bnf] = 2 !!! */ setdown( bnf, bnfext, alpha, embed ) = { local( tmp, deg, leastcoeff, result ); deg = poldegree( lift( alpha ) ); tmp = polcoeff( lift( alpha ), deg )/polcoeff( embed, deg ); leastcoeff = polcoeff( lift( alpha ) - tmp*embed, 0 ); result = Mod( tmp*variable( bnf.pol ) + leastcoeff, bnf.pol ); return( result ); } /* Inverse to setdown */ setup( bnf, bnfext, alpha, embed ) = { local( tmp ); tmp = lift( alpha ); tmp = subst( tmp, variable( bnf.pol ), embed ); tmp = Mod( tmp, bnfext.pol ); return( tmp ); } /* Chinese remainder theorem in the ring of ints of bnf. bnf is supposed to have class number 1. */ china( bnf, X, P ) = { /* darst_ggT( bnf, alpha, beta ) in hnf.gp */ local( i, PE, D, E, myprod, R, tmp ); D = vector( length( P ), i, 0 ); /* Erzeuger von Idealen in P */ R = vector( length( P ), i, 0 ); E = vector( length( P ), i, 0 ); for( i = 1, length( D ), D[i] = nfbasistoalg( bnf, bnfisprincipal( bnf, P[i] )[2] ); ); myprod = 1; for( i = 1, length( D ), myprod = myprod*D[i] ); for( j = 1, length( D ), E[j] = myprod / D[j]; ); for( j = 1, length( D ), tmp = darst_ggT( bnf, D[j], E[j] ); R[j] = E[j]*tmp[2]; ); tmp = 0; for( j = 1, length( R ), tmp = tmp + X[j]*R[j]; ); return( tmp ); } VectorModHNF(v, H)= { local( m, n, i, j, r, q, k ); m = matsize(H)[1]; n = matsize(H)[2]; j = n; i = m; while (j > 0, while (H[i, j] == 0, i--); r = v[i] % H[i, j]; q = (v[i] - r) / H[i, j]; for (k=1, i, v[k] = v[k] - q*H[k,j]); j--; ); return( v ); } /* Inverse of the integral element a mod the ideal Id, where K is supposed to have class number 1. */ InverseModId(K, a, Id)= { local( p, xyd ); p = nfbasistoalg(K, bnfisprincipal(K, Id)[2]); xyd = darst_ggT(K, a, p); if (xyd[3] != 1, print("ModInverse: a nicht teilerfremd zu Id!"); return 0); return( xyd[1] ); } /* Addition of units given in the form [[x1, ..., xm], d], where x1, ..., xm are the coefficients with respect to a fixed system of fundamental S-units. d \in Z is a denominator. */ add( alpha, beta ) = { local( i, tmp, ggT ); if( length( alpha[1] ) != length( beta[1] ), print(" Die Einheiten sind incompatibel " ); return( 0 ); ); for( i = 1, length( alpha[1] ), alpha[1][i] = alpha[1][i]*beta[2] + beta[1][i]*alpha[2]; ); alpha[2] = alpha[2]*beta[2]; ggT = alpha[2]; for( i = 1, length( alpha[1] ), ggT = gcd( alpha[1][i], ggT ); ); for( i = 1, length( alpha[1] ), alpha[1][i] = alpha[1][i]/ggT; ); alpha[2] = alpha[2]/ggT; return( alpha ); } /* Multiplikation of a unit in basis representation with a scalar in Z */ scalmul( alpha, unit ) = { local( i, ggT ); unit[1] = alpha*unit[1]; ggT = unit[2]; for( i = 1, length( unit[1] ), ggT = gcd( unit[1][i], ggT ); ); unit[1] = (1/ggT)*unit[1]; unit[2] = (1/ggT)*unit[2]; return( unit ); } /* Berechnung von Galoisbahn einer Einheit. Liefert eine Matrix, in der aut(unit) als Zeilenvektor der Z-Koordinaten auf Basis- einheiten repraesentiert wird. unit wird als Tensor mit Q ueber- geben. */ unitspan( unit, D, ordG ) = { local( U, i, j, tmp ); U = matrix( ordG - 1, length( unit[1] ), i, j, 0 ); tmp = mattranspose( unit[1] ); for( i = 1, matsize( U )[1], if ( i == 1, for ( j = 1, matsize( U )[2], U[i, j] = tmp[j]; ), tmp = D*tmp; for( j = 1, matsize( U )[2], U[i, j] = tmp[j]; ); ); ); return( U ); } /* Teste, ob eine Einheit trivial ist ( dann ist das Ergebnis 1, sonst 0 ). unit wird als Tensor mit Q uebergeben. */ testnull( unit ) = { local( i, flag ); flag = 1; for( i = 1, length( unit[1]), if( unit[1][i] != 0, flag = 0; return( flag ); ); ); return( flag ); } /* A ist gegeben in HNF. ACHTUNG: RANG( A ) = 2 ! Berechnet werden komplementaere zu diejenigen Zeilen, deren Aufspann zu l teilerfremd ist und minimalen Index besitzt */ findkomp( A, l ) = { local( i, j, k, tmp, puffer, flag, len, ind, n, m ); n = matsize( A )[1]; m = matsize( A )[2]; ind = vector( n, i, i ); len = n!/( m!*( n - m )! ); puffer = vector( len, i ); k = 1; for( i = 1, matsize( A )[1] - 1, for( j = i + 1, matsize( A )[1], tmp = concat( vecextract( A~, [i] ), vecextract( A~, [j] ) ); puffer[k] = [ matdet( tmp ), [i, j] ]; k = k + 1; ); ); k = 1; while( puffer[k][1]%l == 0 && k < len, k = k + 1; ); if( k == len && puffer[k][1]%l == 0, print( " Index ist nicht teilerfremd zu l " ); return( 0 ); ); tmp = abs( puffer[k][1] ); flag = k; i = k; while( i <= len, if( abs( puffer[i][1] ) < tmp && ( puffer[i][1] )%l != 0, tmp = abs( puffer[i][1] ); flag = i; ); i = i + 1; ); return( complement( ind, puffer[flag][2] ) ); } /************************************************************************************************* ********** FUNCTIONS RELATED TO HNF COMPOTATIONS IN PRINCIPAL IDEAL DOMAINS ********************* **************************************************************************************************/ /* ( x -> alpha*x ) bzgl. der Ganzheitsbasis beschreibende Matrix, alpha ist durch polmod gegeben */ mult_mit( bnf, alpha ) = { local( M, i, j, n, tmp, intbasis ); intbasis = bnf.zk; M = matrix( n = poldegree( bnf.pol ), n, i, j ); for( i = 1, n, tmp = nfalgtobasis( bnf, alpha*intbasis[i] ); for( j = 1, n, M[j, i] = tmp[j]; ); ); return( M ); } /* Erweiterte Matrix ( A, B ) bei gegebenen Matrizen A, B */ ext_matrix( A, B ) = { local( M, i, j, n, m, sizeA, sizeB ); sizeA = matsize( A ); sizeB = matsize( B ); if( !( sizeA[1] == sizeB[1] ), print( " Die Matrizen haben unpassende Rahmen " ); return; ); M = matrix( n = sizeA[1], m = sizeA[2] + sizeB[2], i, j ); for( i = 1, n, for( j = 1, m, if( j <= sizeA[2], M[i, j] = A[i, j], M[i, j] = B[i, j - sizeA[2]] ); ); ); return( M ); } /* Matrizenmultiplikation A*B */ matmul( A, B ) = { local( C, i, j, k, n, m, sizeA, sizeB, tmp ); sizeA = matsize( A ); sizeB = matsize( B ); if( !( sizeA[2] == sizeB[1] ), print( " Die Matrizen haben unpassende Rahmen " ); return( 0 ); ); C = matrix( n = sizeA[1], m = sizeB[2], i, j ); for( i = 1, n, for( j = 1, m, tmp = 0; for( k = 1, sizeA[2], tmp = tmp + A[i, k]*B[k, j] ); C[i, j] = tmp; ); ); return( C ); } /* Matrix*Vector */ matvecmul( A, x ) = { local( y, i, j, n, sizeA, tmp ); sizeA = matsize( A ); if( !( sizeA[2] == length( x ) ), print( " Die Matrizen haben unpassende Rahmen " ); return( 0 ); ); y = vector( n = sizeA[2], i ); for( i = 1, sizeA[1], tmp = 0; for( j = 1, n, tmp = tmp + A[i, j]*x[j] ); y[i] = tmp; ); return( mattranspose( y ) ); } /* Multiplikatore zur Darstellung von ggT( x, y ) */ darst_ggT( bnf, alpha, beta ) = { local( i, tmp, A, B, M, X, Y, Z, xi1, xi2, d ); A = mult_mit( bnf, alpha ); B = mult_mit( bnf, beta ); /*return([A,B]);*/ M = mathnf( ext_matrix( A, B ), 1 ); /*return(M);*/ d = bnfisprincipal( bnf, M[1] )[2]; /*return([M, d]);*/ X = matsolve( M[1], d ); tmp = length( X ); Z = vector( 2*tmp, i, 0 ); for( i = tmp + 1, length( Z ), Z[i] = X[i - tmp] ); Z = matvecmul( M[2], Z ); /*Z=M[2]*Z~;*/ /*return( Z );*/ Y = vector( tmp, i ); for( i = 1, tmp, X[i] = Z[i] ); for( i = tmp + 1, length( Z ), Y[i - tmp] = Z[i] ); xi1 = nfbasistoalg( bnf, X ); xi2 = nfbasistoalg( bnf, mattranspose( Y ) ); d = nfbasistoalg( bnf, d ); return( [ xi1, xi2, d ] ); } /* Berechnung eines Repraesentanten modulo Ideal in dem Ring der ganzalgebraischen Zahlen */ rest_mod( bnf, alpha, module ) = { local( i, j, n, m, b, q, f, A ); A = idealhnf( bnf, module ); b = nfalgtobasis( bnf, alpha ); n = matsize( A )[1]; m = matsize( A )[2]; i = n + 1; while( !( i == 1 ), i = i - 1; f = A[i, i]; if( f < 0, for( j = 1, n, A[j, i] = -A[j, i]; ); ); if( !( f == 0 ), q = floor( b[i]/f ); for( j = 1, n, b[j] = b[j] - q*A[j, i]; ); ); ); b = nfbasistoalg( bnf, b ); return( b ); } /* Berechnung der Hermiteschen Normalform einer Matrix A, mit ganzalgebraischen Zahlen aus einem HIR als Eintraegen, Ergebnis ist [ H, U ], wobei H die HNF(A) ist und A*U = H gilt */ mathnf_allg( bnf, A ) = { local( e, i, j, k, l, b, q, n, m, ggt ); local( u, v, d, tmp1, tmp2, U, H, B, B1 ); m = matsize( A )[1]; n = matsize( A )[2]; U = matrix( n, n, i, j, 0 ); B = vector( m, i, 0 ); B1 = vector( n, i, 0 ); for( i = 1, n, U[i, i] = 1 ); i = m + 1; k = n + 1; j = k; if( m <= n, l = 1, l = m - n + 1; ); while( !( i == l ), i = i - 1; k = k - 1; j = k; while( !( j == 1 ), j = j - 1; if( !( A[i,j] == 0 ), /* Euklidean Step */ ggt = darst_ggT( bnf, A[i, k], A[i, j] ); /*if( i == i, print( "j =", j,"i=", i,"k=", k,"ggT = ", ggt ); );*/ u = ggt[1]; v = ggt[2]; d = ggt[3]; tmp1 = A[i, k]; tmp2 = A[i, j]; for( e = 1, m, B[e] = u*A[e, k] + v*A[e, j] ); for( e = 1, m, A[e, j] = tmp1/d*A[e, j] - tmp2/d*A[e, k] ); for( e = 1, m, A[e,k] = B[e] ); for( e = 1, n, B1[e] = u*U[e, k] + v*U[e, j] ); for( e = 1, n, U[e, j] = tmp1/d*U[e, j] - tmp2/d*U[e, k] ); for( e = 1, n, U[e,k] = B1[e] ); ); ); /*return( [A, U] );*/ b = A[i, k]; /* Final Reduction */ if( b == 0, k = k + 1, for( j = k + 1, n, tmp1 = A[i, j]; tmp2 = rest_mod( bnf, tmp1, d ); tmp1 = tmp1 - tmp2; q = tmp1/b; for( e = 1, m, A[e, j] = A[e, j] - q*A[e, k]; ); for( e = 1, n, U[e, j] = U[e, j] - q*U[e, k]; ); ); ); ); return( [A, U] ); } /* Apply iso to the algebraic number a. */ {Sigma(a, minpoly, iso)= Mod(subst(lift(a), x, iso), minpoly); } /* Apply iso^n to the algebraic number a. */ {Sigmapow(a, minpoly, iso, n, i, z) = z = a; for (i=1, n, z = Sigma(z, minpoly, iso)); return(z); } {VecGcd(v, n, d, i, j) = n = length(v); if ( n == 1, return(v[1]) ); d = v[1]; for (i=1, n, d = gcd(d, v[i])); return(d); } /* Fuer zwei Mengen M, N, gegeben als Vektoren berechnet die mengentheoretische Differenz M\N */ complement( M, N ) = { local( puffer, i, j, mlen, nlen, tmp, flag ); puffer = List( [] ); mlen = length( M ); nlen = length( N ); for( j = 1, mlen, flag = 1; tmp = M[j]; i = 1; while( flag != 0 && i <= nlen, if( tmp == N[i], flag = 0, ); i = i + 1; ); if( flag == 1, puffer = concat( puffer, List( [tmp] ) ); ); ); return( puffer ); } /* Potenzen einer Zahl modulo Ideal */ power( bnf, alpha, n, module ) = { local( pow, len ); if( n == 0, return ( 1 ); ); pow = alpha; len = length( binary( n ) ) - 1; while( !( len == 0 ), len = len - 1; pow = pow*pow; pow = rest_mod( bnf, pow, module ); if( bittest( n, len ), pow = pow*alpha; pow = rest_mod( bnf, pow, module ); ); ); return( pow ); } /* Betragsgroesste Fehler beim Runden der Koeffizienten eines Polynoms */ absdiff( pol ) = { local( i, roundpol, coeff, roundcoeff, diff, tmp, deg ); deg = poldegree( pol ); coeff = vector( deg + 1, i, 0 ); roundcoeff = vector( deg + 1, i, 0 ); roundpol = round( real( pol ) ); for( i = 0, deg, coeff[i + 1] = polcoeff( pol, i ); roundcoeff[i + 1 ] = polcoeff( roundpol, i ); ); diff = abs( coeff[1] - roundcoeff[1] ); for( i = 2, length( coeff ), tmp = abs( coeff[i] - roundcoeff[i] ); if( tmp > diff, diff = tmp; ); ); return( diff ); }