diff --git a/src/Map/NfOrd.jl b/src/Map/NfOrd.jl index 5d01396592..021f35fbc1 100644 --- a/src/Map/NfOrd.jl +++ b/src/Map/NfOrd.jl @@ -1087,6 +1087,7 @@ function image(mF::NfToGFMor_easy, a::FacElem{AbsSimpleNumFieldElem, AbsSimpleNu p = mF.defining_pol q = one(Fq) t = mF.t + quo = gcd(quo, order(Fq)-1) for (k, v) = a.fac vv = v if quo != 0 diff --git a/src/NumField/NfAbs.jl b/src/NumField/NfAbs.jl index 4aa2f43f60..93f456e9b3 100644 --- a/src/NumField/NfAbs.jl +++ b/src/NumField/NfAbs.jl @@ -21,3 +21,4 @@ include("NfAbs/MPolyFactor.jl") include("NfAbs/MPolyAbsFact.jl") include("NfAbs/ConjugatesNS.jl") include("NfAbs/Cyclotomic.jl") +include("NfAbs/MultDep.jl") diff --git a/examples/MultDep.jl b/src/NumField/NfAbs/MultDep.jl similarity index 53% rename from examples/MultDep.jl rename to src/NumField/NfAbs/MultDep.jl index dfe934b2bf..2c30422420 100644 --- a/examples/MultDep.jl +++ b/src/NumField/NfAbs/MultDep.jl @@ -3,11 +3,29 @@ module MultDep using Hecke import Base.* -function multiplicative_group_mod_units_fac_elem(A::Vector{AbsSimpleNumFieldElem}; use_max_ord::Bool = false) +""" +Given A[i] elements in K, find matrices I and U s.th. +A^I is a basis for + /Units < K^*/Units +actually a sub-group of the S-unit group for a coprime base +generated from the elements in A. + +A^U is a generating set for the relations: + < Units +and every u s.th. A^u is a unit is in the span of U + +The coprime base is returned as well, the columns of I correspond to the +ordeing of the base. +""" +function syzygies_sunits_mod_units(A::Vector{AbsSimpleNumFieldElem}; use_ge::Bool = false, max_ord::Union{Nothing, AbsSimpleNumFieldOrder}=nothing) k = parent(A[1]) @assert all(i->parent(i) === k, A) - if use_max_ord - zk = maximal_order(k) + if !use_ge + if max_ord === nothing + zk = maximal_order(k) + else + zk = max_ord + end cp = coprime_base(A, zk) else cp = coprime_base(A) @@ -26,10 +44,20 @@ function multiplicative_group_mod_units_fac_elem(A::Vector{AbsSimpleNumFieldElem end push!(M, sparse_row(FlintZZ, T)) end + h, t = Hecke.hnf_with_transform(matrix(M)) + h = h[1:rank(h), :] + return h, t[nrows(h)+1:end, :], cp + # THINK! do we want or not... + # - M is naturally sparse, hence it makes sense + # - for this application we need all units, hence the complete + # kernel - which is huge and dense... + # - cp seems to not be used for anything afterwards. + # It will be when we actually create the group, in the DiscLog h, t = Hecke.hnf_kannan_bachem(M, Val(true), truncate = true) return h, t, cp end +#= function units(h::SMat, t, b::Vector{AbsSimpleNumFieldElem}) u = FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField}[] for i=nrows(h)+1:length(b) @@ -42,51 +70,16 @@ function units(h::SMat, t, b::Vector{AbsSimpleNumFieldElem}) end return u end +=# -function unit_group_mod_torsion_fac_elem(O::AbsNumFieldOrder, u::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField}}) - U = Hecke._unit_group_init(O) - s = signature(O) - r = s[1] + s[2] - 1 - for y = u - is_tors, p = is_torsion_unit(y, false, U.tors_prec) - U.tors_prec = max(p, U.tors_prec) - if is_tors - continue - end - Hecke._add_unit(U, y) - if length(U.units) >= r - break - end - end - if length(U.units) < r - # maybe use pAdic stuff here... - error("not complete yet") - end - - U.full_rank = true - - U.units = Hecke.reduce(U.units, U.tors_prec) - - for y = u - is_tors, p = is_torsion_unit(y, false, U.tors_prec) - U.tors_prec = max(p, U.tors_prec) - if is_tors - continue - end - x = Hecke.reduce_mod_units([y], U)[1] - is_tors, p = is_torsion_unit(x, false, U.tors_prec) - U.tors_prec = max(p, U.tors_prec) - if is_tors - continue - end - Hecke._add_dependent_unit(U, x) - end - return U -end - +#= For Ge's PhD + - smallest common overorder (non-coprime case) + - ideals in non-maximal orders with checking if the operation + worked - if not automatically redoing with larger order +=# function *(O1::AbsNumFieldOrder, O2::AbsNumFieldOrder) - k = nf(O1) - @assert k === nf(O2) + k = number_field(O1) + @assert k === number_field(O2) b1 = basis(O1, k) n = length(b1) b2 = basis(O2, k) @@ -169,7 +162,11 @@ function divexact(a::GeIdeal, b::GeIdeal; check::Bool=true) end Hecke.norm(a::GeIdeal) = norm(a.a) - +#XXX: not sure if those 2 should get "exported", Ge does not seem to be +# overly fast.... +#TODO: do an integer coprime base first, then refine. Possibly also do +# a p-maximality for all p from the integer coprime base. +#TODO: find examples where Ge is useful function coprime_base(A::Vector{AbsSimpleNumFieldElem}) c = Vector{GeIdeal}() for a = A @@ -177,7 +174,7 @@ function coprime_base(A::Vector{AbsSimpleNumFieldElem}) isone(n) || push!(c, n) isone(d) || push!(c, d) end - return coprime_base(c) + return Hecke.AbstractAlgebra.coprime_base_steel(c) end function coprime_base(A::Vector{AbsSimpleNumFieldElem}, O::AbsNumFieldOrder) @@ -195,8 +192,29 @@ function valuation(a::AbsSimpleNumFieldElem, p::GeIdeal) return valuation(a, p.a) end +""" +reduce the rows of v modulo the lattice spanned by the rows of M. +M should be LLL reduced. +""" +function size_reduce(M::ZZMatrix, v::ZZMatrix) + s = gram_schmidt_orthogonalisation(QQ.(transpose(M))) + d = diagonal(transpose(s)*s) + for i=1:nrows(v) + for j=nrows(M):-1:1 + x = ZZ(round((v[i:i, :]* s[:, j:j]/d[j])[1,1])) + v[i:i, :] -= x*M[j:j, :] + end + end + return v +end -function mult_syzygies_units(A::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField}}) +""" +A a vector of units in fac-elem form. Find matrices U and V s.th. +A^U is a basis for /Tor +and +A^V is a generating system for the relations of A in Units/Tor +""" +function syzygies_units_mod_tor(A::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField}}) p = next_prime(100) K = base_ring(parent(A[1])) m = maximum(degree, keys(factor(GF(p), K.pol).fac)) @@ -211,16 +229,22 @@ function mult_syzygies_units(A::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleN r1, r2 = signature(K) r = r1+r2 -1 n = degree(K) - C = qAdicConj(K, p) - la = matrix(conjugates_log(A[1], C, prec, all = false, flat = true))' + C = Hecke.qAdicConj(K, p) + la = transpose(matrix(conjugates_log(A[1], C, prec, all = false, flat = true))) lu = zero_matrix(base_ring(la), 0, n) uu = [] - for a = A + indep = Int[] + dep = Int[] + for ia = 1:length(A) + a = A[ia] while true - @vtime :qAdic 1 la = matrix(conjugates_log(a, C, prec, all = false, flat = true))' + @vtime :qAdic 1 la = transpose(matrix(conjugates_log(a, C, prec, all = false, flat = true))) if iszero(la) @vtime :qAdic 1 @hassert :qAdic 1 verify_gamma([a], [ZZRingElem(1)], ZZRingElem(p)^prec) @vprint :qAdic 1 println("torsion found") + gamma = vcat([ZZ(0) for i=1:r+length(uu)], [ZZ(1)]) + push!(uu, (a, gamma)) + push!(dep, ia) break end lv = vcat(lu, la) @@ -237,18 +261,20 @@ function mult_syzygies_units(A::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleN @vtime :qAdic 1 k = kernel(lv, side = :left) @assert nrows(k) < 2 if nrows(k) == 0 - println("new ") + @vprint :qAdic 1 "new independent unit found\n" push!(u, a) + push!(indep, ia) lu = vcat(lu, la) @assert length(u) <= r else # length == 1 extend the module + @vprint :qAdic 1 "dependent unit found, looking for relation\n" s = QQFieldElem[] for x in k[1, :] @vtime :qAdic 1 y = lift_reco(FlintQQ, x, reco = true) if y === nothing prec *= 2 @vprint :qAdic 1 "increase prec to ", prec - lu = matrix([conjugates_log(x, C, prec, all = false, flat = true) for x = u])' + lu = transpose(matrix([conjugates_log(x, C, prec, all = false, flat = true) for x = u])) break end push!(s, y) @@ -259,7 +285,7 @@ function mult_syzygies_units(A::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleN d = reduce(lcm, map(denominator, s)) gamma = ZZRingElem[FlintZZ(x*d)::ZZRingElem for x = s] @assert reduce(gcd, gamma) == 1 # should be a primitive relation - @time if !verify_gamma(push!(copy(u), a), gamma, ZZRingElem(p)^prec) + if !verify_gamma(push!(copy(u), a), gamma, ZZRingElem(p)^prec) prec *= 2 @vprint :qAdic 1 "increase prec to ", prec lu = matrix([conjugates_log(x, C, prec, all = false, flat = true) for x = u])' @@ -268,6 +294,7 @@ function mult_syzygies_units(A::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleN @assert length(gamma) == length(u)+1 gamma = vcat(gamma[1:length(u)], [0 for i=length(u)+1:r+length(uu)], [gamma[end]]) push!(uu, (a, gamma)) + push!(dep, ia) end break end @@ -277,14 +304,19 @@ function mult_syzygies_units(A::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleN has rank s and r_i in Z^n be such that prod u_i^r_i = 1 (OK, sum of the logs is zero) - rank = s as well - so the r_i form a Q-basis for the relations. + rank = s as well #wrong: + #Input: u^2, u^3, u^7 in rank r = 2 + #U = [-3 0 2 0; -7 0 0 2] + the r_i form a Q-generating basis for the relations. + They are indep as they are growing in length (new columns per new element, non + are torsion, so the final non-zero entries are moving right.) Essentially, the gamma of above are the r_i Let [H|0] = [r_i|i]*T be the hnf with trafo, so T in Gl(n, Z) Then = <[u_i|i] T> [r_i|i] * [u_i|i]^t = 0 (by construction) [r_i|i] T inv(T) [u[i] | i] = 0 + Careful! H may not have full rank, hence the shape is different [H | 0] [v_i | i] = 0 so, since H is triangular(!!) v_1, ... v_n-s = 0 and = @@ -292,35 +324,48 @@ function mult_syzygies_units(A::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleN for the case of n=s+1 this is mostly the "normal" construction. Note: as a side, the relations do not have to be primitive. If they are, (and n=s+1), then H = 1 + We deal with that by saturation =# for i=1:length(uu)-1 append!(uu[i][2], zeros(FlintZZ, length(uu[end][2])-length(uu[i][2]))) end - if length(uu) == 0 - U = matrix(FlintZZ, length(uu), length(uu[end][2]), reduce(vcat, [x[2] for x = uu])) + if length(uu) == 0 #all torsion + return [], A else U = matrix(FlintZZ, length(uu), length(uu[end][2]), reduce(vcat, [x[2] for x = uu])) + U = hcat(U[:, 1:length(u)], U[:, r+1:ncols(U)]) end - _, U = hnf_with_transform(U') - if false - U = inv(U) - V = sub(U, 1:rows(U), 1:cols(U)-length(u)) - U = sub(U, 1:rows(U), cols(U)-length(u)+1:cols(U)) - #U can be reduced modulo V... - Z = zero_matrix(FlintZZ, cols(V), cols(U)) - I = identity_matrix(FlintZZ, cols(U)) * p^(2*prec) - k = base_ring(A[1]) - A = [ Z V'; I U'] - l = lll(A) - U = sub(l, cols(V)+1:rows(l), cols(U)+1:cols(l)) - U = lll(U) - else - U = lll(U') + + U = saturate(U) + + _, U = hnf_with_transform(transpose(U)) + + k = base_ring(A[1]) + + U = inv(U) + V = sub(U, 1:nrows(U), 1:ncols(U)-length(u)) #the torsion part + U = sub(U, 1:nrows(U), ncols(U)-length(u)+1:ncols(U)) + #U can be reduced modulo V... + V = lll(transpose(V)) + U = size_reduce(V, transpose(U)) + U = lll(U) + + #so basis is A[indep] cat A[dep] ^U + #rels: A[tor], .. * V + nt = zero_matrix(ZZ, length(A), length(A)) + for i=1:length(indep) + nt[indep[i], i] = 1 + end + for i=1:length(dep) + nt[dep[i], i+length(indep)] = 1 end - return Hecke._transform(vcat(u, FacElem{AbsSimpleNumFieldElem,AbsSimpleNumField}[FacElem(k(1)) for i=length(u)+1:r], [x[1] for x = uu]), U') + @assert matrix([collect(1:length(A))]) * nt == matrix([vcat(indep, dep)]) + rel = nt*transpose(V) + return nt*transpose(U), rel end + function verify_gamma(a::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField}}, g::Vector{ZZRingElem}, v::ZZRingElem) #knowing that sum g[i] log(a[i]) == 0 mod v, prove that prod a[i]^g[i] is #torsion @@ -338,7 +383,7 @@ function verify_gamma(a::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField # and, see the bottom, \|Log()\|_2^2 >= 1/4 arcosh((B-2)/2)^2 B = ArbField(nbits(v)*2)(v)^2 B = 1/2 *acosh((B-2)/2)^2 - p = Hecke.upper_bound(log(B)/log(parent(B)(2)), ZZRingElem) + p = Hecke.upper_bound(ZZRingElem, log(B)/log(parent(B)(2))) @vprint :qAdic 1 "using", p, nbits(v)*2 b = conjugates_arb_log(t, max(-Int(div(p, 2)), 2)) # @show B , sum(x*x for x = b), is_torsion_unit(t)[1] @@ -346,9 +391,25 @@ function verify_gamma(a::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField return B > sum(x*x for x = b) end +""" +A padic number a is internally written as + p^v*u +for a unit u given in some precision. +This returns u, v and the precision +""" +function canonical_split(a::PadicFieldElem) + u = ZZ() + ccall((:fmpz_set, Hecke.libflint), Nothing, (Ref{ZZRingElem}, Ref{Int}), u, a.u) + return u, a.v, a.N +end + +#TODO: different name ... function lift_reco(::QQField, a::PadicFieldElem; reco::Bool = false) + if iszero(a) + return QQ(0) + end if reco - u, v, N = getUnit(a) + u, v, N = canonical_split(a) R = parent(a) fl, c, d = rational_reconstruction(u, prime(R, N-v)) !fl && return nothing @@ -360,10 +421,84 @@ function lift_reco(::QQField, a::PadicFieldElem; reco::Bool = false) return x*prime(R, v) end else - return lift(FlintQQ, a) + return lift(QQ, a) end end -Hecke.number_of_rows(A::Matrix{T}) where {T} = size(A)[1] -Hecke.number_of_columns(A::Matrix{T}) where {T} = size(A)[2] +""" +Given torsion units A[i] in factored form, find a generator +for the group as well as a basis for the relations +""" +function syzygies_tor(A::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField}}) + K = base_ring(parent(A[1])) + bnd = Int(Hecke._torsion_group_order_divisor(K)) + #bnd is a multiple of the torsion in the field + #so any residue field where bnd divides q-1 can be used to find the + #relations, Tor is a subgroup of k^* in this case. + O = any_order(K) + #XXX: bad, but prime ideals in non-max orders seem a bit dodgy + # in particular extend_easy gives garbage + #= + k, a = wildanger_field(5,13); + O = any_order(k); + P = prime_ideals_over(O, 71)[1] + F, mF = residue_field(O, P) + mF(O(gen(k))) # => 25 + mmF = Hecke.extend_easy(mF, k); + mmF(gen(k)) # => 70 + degree(P) # => 0 + P.prim_elem # => 21 but should be a root... + =# + + O = maximal_order(K) + for p = PrimesSet(bnd, -1, bnd, 1) + if discriminant(O) % p == 0 + continue + end + P = prime_ideals_over(O, p)[1] + if degree(P) > 1 + continue + end + F, mF = Hecke.ResidueFieldSmallDegree1(O, P) + mF = Hecke.extend_easy(mF, K) + a = try + map(mF, A) + catch e + if isa(e, Hecke.BadPrime) + #if some element in the base of the fac elems is not a unit mod P + #use the next prime... + continue + end + rethrow(e) + end + U, mU = unit_group(F) + b = map(pseudo_inv(mU), a) + B = free_abelian_group(length(A)) + h = hom(B, U, b) + k, mk = kernel(h) + i, mi = image(h) + @assert ngens(i) == 1 + return preimage(h, mi(i[1])).coeff, vcat([mk(x).coeff for x = gens(k)]...) + end +end + +""" +Given non-zero elements A[i] in K, find a basis for the relations, returned + as a matrix. +""" +function syzygies(A::Vector{AbsSimpleNumFieldElem}; use_ge::Bool = false, max_ord::Union{Nothing, AbsSimpleNumFieldOrder} = nothing) + _, t, _ = syzygies_sunits_mod_units(A; use_ge, max_ord) + u = [FacElem(A, t[i, :]) for i = 1:nrows(t)] + _, tt = syzygies_units_mod_tor(u) + u = Hecke._transform(u, tt) + _, ttt = syzygies_tor(u) + return ttt*transpose(tt)*t +end + +export syzygies + +end + +using .MultDep + diff --git a/src/Sparse/UpperTriangular.jl b/src/Sparse/UpperTriangular.jl index 830ff289d4..889c65d955 100644 --- a/src/Sparse/UpperTriangular.jl +++ b/src/Sparse/UpperTriangular.jl @@ -92,7 +92,7 @@ function elementary_divisors(A::SMat{ZZRingElem}) e[i] *= p^v[i] end end - for i = length(e):min(nrows(A), ncols(A)) + for i = length(e)+1:min(nrows(A), ncols(A)) push!(e, ZZRingElem(0)) end return e diff --git a/test/NumField/NfAbs/NfAbs.jl b/test/NumField/NfAbs/NfAbs.jl index 8b3460f127..1e32a3c5dc 100644 --- a/test/NumField/NfAbs/NfAbs.jl +++ b/test/NumField/NfAbs/NfAbs.jl @@ -85,3 +85,17 @@ end @test !is_normal(K) @test length(automorphism_list(K)) == 2 end + + +@testset "NumField/NfAbs/MultDep" begin + k, a = wildanger_field(5,13); + zk = lll(maximal_order(k)) + class_group(zk) + h = zk.__attrs[:ClassGrpCtx] + r = vcat(h.R_gen, h.R_rel); + r = [x for x = r if isa(x, AbsSimpleNumFieldElem)] + q = Hecke.syzygies(r) + @test all(isone, evaluate(FacElem(r, q[i, :])) for i=1:nrows(q)) +end + +