<commit_before>################################################################################
#
#       GrpAb/GrpAbFinGen.jl : Finitely generated abelian groups
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: <NAME>, <NAME>
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 <NAME>, <NAME>
#
################################################################################

export abelian_group, free_abelian_group, issnf, ngens, nrels, rels, snf, isfinite,
       isinfinite, rank, order, exponent, istrivial, isisomorphic,
       direct_product, istorsion, torsion_subgroup, sub, quo, iscyclic,
       psylow_subgroup, issubgroup, abelian_groups, flat, tensor_product,
       dual, chain_complex, isexact, homology, free_resolution, obj, map,
       primary_part, isfree

import Base.+, Nemo.snf, Nemo.parent, Base.rand, Nemo.issnf

################################################################################
#
#  Functions for some abstract interfaces
#
################################################################################

elem_type(G::GrpAbFinGen) = GrpAbFinGenElem

elem_type(::Type{GrpAbFinGen}) = GrpAbFinGenElem

parent_type(::Type{GrpAbFinGenElem}) = GrpAbFinGen

isabelian(::GrpAbFinGen) = true

##############################################################################
#
#  Constructors
#
##############################################################################

@doc Markdown.doc"""
    abelian_group(::Type{T} = GrpAbFinGen, M::fmpz_mat) -> GrpAbFinGen

Creates the abelian group with relation matrix `M`. That is, the group will
have `ncols(M)` generators and each row of `M` describes one relation.
"""
abelian_group(M::fmpz_mat; name::String = "") = abelian_group(GrpAbFinGen, M, name=name)

function abelian_group(::Type{GrpAbFinGen}, M::fmpz_mat; name::String = "")
   if issnf(M) && nrows(M) > 0  && ncols(M) > 0 && !isone(M[1, 1])
    N = fmpz[M[i, i] for i = 1:min(nrows(M), ncols(M))]
    if ncols(M) > nrows(M)
      N = vcat(N, fmpz[0 for i = 1:ncols(M)-nrows(M)])
    end
    G = GrpAbFinGen(N)
  else
    G = GrpAbFinGen(M)
  end
  name == "" || set_name!(G, name)
  return G
end

@doc Markdown.doc"""
    abelian_group(::Type{T} = GrpAbFinGen, M::AbstractMatrix{<:IntegerUnion})

Creates the abelian group with relation matrix `M`. That is, the group will
have `ncols(M)` generators and each row of `M` describes one relation.
"""
function abelian_group(M::AbstractMatrix{<:IntegerUnion}; name::String = "")
  return abelian_group(GrpAbFinGen, M, name=name)
end

function abelian_group(::Type{GrpAbFinGen}, M::AbstractMatrix{<:IntegerUnion}; name::String = "")
  return abelian_group(matrix(FlintZZ, M), name=name)
end

function _issnf(N::Vector{T}) where T <: IntegerUnion
  if isone(length(N)) && isone(N[1])
    return false
  end
  for i = 1:length(N)-1
    if isone(abs(N[i]))
      return false
    end
    if iszero(N[i])
      if !iszero(N[i+1])
        return false
      end
    elseif !iszero(mod(N[i+1], N[i]))
      return false
    end
  end
  return true
end

@doc Markdown.doc"""
    abelian_group(::Type{T} = GrpAbFinGen, M::AbstractVector{<:IntegerUnion}) -> GrpAbFinGen
    abelian_group(::Type{T} = GrpAbFinGen, M::IntegerUnion...) -> GrpAbFinGen

Creates the direct product of the cyclic groups $\mathbf{Z}/m_i$,
where $m_i$ is the $i$th entry of `M`.
"""
function abelian_group(M::AbstractVector{<:IntegerUnion}; name::String = "")
  return abelian_group(GrpAbFinGen, M, name=name)
end

function abelian_group(::Type{GrpAbFinGen}, M::AbstractVector{<:IntegerUnion}; name::String = "")
  if _issnf(M)
    G = GrpAbFinGen(M)
  else
	N = zero_matrix(FlintZZ, length(M), length(M))
    for i = 1:length(M)
      N[i,i] = M[i]
    end
    G = GrpAbFinGen(N)
  end
  if !isempty(M)
    res = lcm(M)
    if !iszero(res)
      G.exponent = res
    end
  end
  name == "" || set_name!(G, name)
  return G
end

function abelian_group(M::IntegerUnion...; name::String = "")
  return abelian_group(collect(M), name=name)
end

@doc Markdown.doc"""
    free_abelian_group(::Type{T} = GrpAbFinGen, n::Int) -> GrpAbFinGen

Creates the free abelian group of rank `n`.
"""
free_abelian_group(n::Int) = free_abelian_group(GrpAbFinGen, n)

function free_abelian_group(::Type{GrpAbFinGen}, n::Int)
  return abelian_group(GrpAbFinGen, zeros(Int, n))
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, A::GrpAbFinGen)
  @show_name(io, A)
  @show_special(io, A)

  if issnf(A)
    show_snf(io, A)
  else
    show_gen(io, A)
  end
end

function show_hom(io::IO, G)#::GrpAbFinGen)
  D = get_special(G, :hom)
  D === nothing && error("only for hom")
  print(io, "hom of ")
  print(IOContext(io, :compact => true), D)
end

function show_direct_product(io::IO, G)#::GrpAbFinGen)
  D = get_special(G, :direct_product)
  D === nothing && error("only for direct products")
  print(io, "direct product of ")
  show(IOContext(io, :compact => true), D)
end

function show_tensor_product(io::IO, G)#::GrpAbFinGen)
  D = get_special(G, :tensor_product)
  D === nothing && error("only for tensor products")
  print(io, "tensor product of ")
  show(IOContext(io, :compact => true), D)
end

function show_gen(io::IO, A::GrpAbFinGen)
  print(io, "(General) abelian group with relation matrix\n$(A.rels)")
  if isdefined(A, :snf_map)
    println(io, "\nwith structure of ", domain(A.snf_map))
  end
end

function show_snf(io::IO, A::GrpAbFinGen)
  if get(io, :compact, false)
    print(io, "Abelian group with structure: ")
  else
    print(io, "GrpAb: ")
  end
  show_snf_structure(io, A)
end

function show_snf_structure(io::IO, A::GrpAbFinGen, mul = "x")
  @assert issnf(A)
  len = length(A.snf)
  if len == 0
    print(io, "Z/1")
    return
  end

  if A.snf[1] == 0
    if len == 1
      print(io, "Z")
    else
      print(io, "Z^$(len)")
    end
    return
  end

  i = 1
  while i <= len
    inv = A.snf[i]
    j = 1
    while i + j <= len && inv == A.snf[i + j]
      j += 1
    end
    if iszero(inv)
      print(io, "Z")
      if j > 1
        print(io, "^($j)")
      end
    else
      if j > 1
        print(io, "(Z/$(inv))^$(j)")
      else
        print(io, "Z/$(inv)")
      end
    end
    if i + j - 1 < len
      print(io, " $mul ")
    end
    i += j
  end
end

################################################################################
#
#  Hash function
#
################################################################################

hash_snf(A::GrpAbFinGen, h::UInt) = hash(A.snf, h)

function hash(A::GrpAbFinGen, h::UInt)
  if issnf(A)
    return hash_snf(A, h)
  else
    return hash_snf(snf(A)[1], h)
  end
end

################################################################################
#
#  Field access
#
################################################################################

@doc Markdown.doc"""
    issnf(G::GrpAbFinGen) -> Bool

Returns whether the current relation matrix of the group $G$ is in Smith
normal form.
"""
issnf(A::GrpAbFinGen) = A.issnf

@doc Markdown.doc"""
    ngens(G::GrpAbFinGen) -> Int

Returns the number of generators of $G$ in the current representation.
"""
function ngens(A::GrpAbFinGen)
  if issnf(A)
    return length(A.snf)
  else
    return ncols(A.rels)
  end
end

@doc Markdown.doc"""
    nrels(G::GrpAbFinGen) -> Int

Returns the number of relations of $G$ in the current representation.
"""
function nrels(A::GrpAbFinGen)
  if issnf(A)
    return length(A.snf)
  else
    return nrows(A.rels)
  end
end

@doc Markdown.doc"""
    rels(A::GrpAbFinGen) -> fmpz_mat

Returns the currently used relations of $G$ as a single matrix.
"""
function rels(A::GrpAbFinGen)
  if issnf(A)
    return rels_snf(A)
  else
    return rels_gen(A)
  end
end

function rels_gen(A::GrpAbFinGen)
  return A.rels
end

function rels_snf(A::GrpAbFinGen)
  M = zero_matrix(FlintZZ, ngens(A), ngens(A))
  for i = 1:ngens(A)
    M[i,i] = A.snf[i]
  end
  return M
end

################################################################################
#
#  Hermite normal form
#
################################################################################

function assure_has_hnf(A::GrpAbFinGen)
  if isdefined(A, :hnf)
    return nothing
  end
  if isdefined(A, :exponent) && nrows(A.rels) >= ncols(A.rels)
    A.hnf = hnf_modular_eldiv(A.rels, A.exponent)
  else
    A.hnf = hnf(A.rels)
  end
  return nothing
end

################################################################################
#
#  Smith normal form
#
################################################################################

@doc Markdown.doc"""
    snf(A::GrpAbFinGen) -> GrpAbFinGen, Map

Returns a pair $(G, f)$, where $G$ is an abelian group in canonical Smith
normal form isomorphic to $G$ and an isomorphism $f : G \to A$.
"""
function snf(G::GrpAbFinGen)
  if isdefined(G, :snf_map)
    return domain(G.snf_map)::GrpAbFinGen, G.snf_map::GrpAbFinGenMap
  end

  if issnf(G)
    G.snf_map = GrpAbFinGenMap(G) # identity
    return G, G.snf_map::GrpAbFinGenMap
  end
  if isdefined(G, :exponent)
    if isdefined(G, :hnf)
      S, T = snf_for_groups(G.hnf, G.exponent)
    else
      S, T = snf_for_groups(G.rels, G.exponent)
    end
  else
    S, _, T = snf_with_transform(G.rels, false, true)
  end

  m = min(nrows(S), ncols(S))
  if m > 0 && nrows(S) >= ncols(S)
    e = S[m, m]
    if e > 1
      if fits(Int, e) && isprime(e)
        F = GF(Int(e), cached = false)
        TF = map_entries(F, T)
        iT = lift(inv(TF))
      else
        iT = invmod(T, e)
      end
    else
      iT = inv(T)
    end
  else
    iT = inv(T)
  end
  return _reduce_snf(G, S, T, iT)
end

# For S in SNF with G.rels = U*S*T and Ti = inv(T) this removes
# the ones at the diagonal of S and constructs the homomorphism.
function _reduce_snf(G::GrpAbFinGen, S::fmpz_mat, T::fmpz_mat, Ti::fmpz_mat)

  d = fmpz[S[i,i] for i = 1:min(nrows(S), ncols(S))]

  while length(d) < ngens(G)
    push!(d, 0)
  end

  pos = Int[i for i = 1:length(d) if !isone(d[i])]
  r = Int[i for i = 1:nrows(T)]
  s = fmpz[ d[i] for i in pos]
  TT = sub(T, r, pos)
  TTi = sub(Ti, pos, r)

  H = GrpAbFinGen(s)
  if !isempty(s) && !iszero(s[end])
    H.exponent = s[end]
    G.exponent = s[end]
  end
  mp = hom(H, G, TTi, TT, check = false)
  G.snf_map = mp
  return H, mp::GrpAbFinGenMap
end

################################################################################
#
#  Finiteness
#
################################################################################

@doc Markdown.doc"""
    isfinite(A::GrpAbFinGen) -> Bool

Returns whether $A$ is finite.
"""
isfinite(A::GrpAbFinGen) = issnf(A) ? isfinite_snf(A) : isfinite_gen(A)

isfinite_snf(A::GrpAbFinGen) = length(A.snf) == 0 || !iszero(A.snf[end])

isfinite_gen(A::GrpAbFinGen) = isfinite(snf(A)[1])

@doc Markdown.doc"""
    isinfinite(A::GrpAbFinGen) -> Bool

Returns whether $A$ is infinite.
"""
isinfinite(A::GrpAbFinGen) = !isfinite(A)

################################################################################
#
#  Rank
#
################################################################################

@doc Markdown.doc"""
    rank(A::GrpAbFinGen) -> Int

Returns the rank of $A$, that is, the dimension of the
$\mathbf{Q}$-vectorspace $A \otimes_{\mathbf Z} \mathbf Q$.
"""
rank(A::GrpAbFinGen) = issnf(A) ? rank_snf(A) : rank_gen(A)

rank_snf(A::GrpAbFinGen) = length(findall(x -> iszero(x), A.snf))

rank_gen(A::GrpAbFinGen) = rank(snf(A)[1])

rank(A::GrpAbFinGen, p::Union{Int, fmpz}) = issnf(A) ? rank_snf(A, p) : rank_snf(snf(A)[1], p)

function rank_snf(A::GrpAbFinGen, p::Union{Int, fmpz})
  if isempty(A.snf)
    return 0
  end
  if !iszero(mod(A.snf[end], p))
    return 0
  end
  i = findfirst(x -> iszero(mod(x, p)), A.snf)
  return length(A.snf)-i+1
end


################################################################################
#
#  Order
#
################################################################################

@doc Markdown.doc"""
    order(A::GrpAbFinGen) -> fmpz

Returns the order of $A$. It is assumed that $A$ is finite.
"""
order(A::GrpAbFinGen) = issnf(A) ? order_snf(A) : order_gen(A)

function order_snf(A::GrpAbFinGen)
  isinfinite(A) && error("Group must be finite")
  return prod(A.snf)
end

order_gen(A::GrpAbFinGen) = order(snf(A)[1])

################################################################################
#
#  Exponent
#
################################################################################

@doc Markdown.doc"""
    exponent(A::GrpAbFinGen) -> fmpz

Returns the exponent of $A$. It is assumed that $A$ is finite.
"""
function exponent(A::GrpAbFinGen)
  if issnf(A)
    res = exponent_snf(A)
    if !iszero(res)
      A.exponent = res
    end
    return res
  else
    res = exponent_gen(A)
    if !iszero(res)
      A.exponent = res
    end
    return res
  end
end

function exponent_snf(A::GrpAbFinGen)
  isinfinite(A) && error("Group must be finite")
  ngens(A)==0 && return fmpz(1)
  return A.snf[end]
end

exponent_gen(A::GrpAbFinGen) = exponent(snf(A)[1])

################################################################################
#
#  Trivial
#
################################################################################

@doc Markdown.doc"""
    istrivial(A::GrpAbFinGen) -> Bool

Checks if $A$ is the trivial group.
"""
istrivial(A::GrpAbFinGen) = isfinite(A) && isone(order(A))

################################################################################
#
#  Isomorphism
#
################################################################################

@doc Markdown.doc"""
    isisomorphic(G::GrpAbFinGen, H::GrpAbFinGen) -> Bool

Checks if $G$ and $H$ are isomorphic.
"""
function isisomorphic(G::GrpAbFinGen, H::GrpAbFinGen)
  b = filter(x -> x != 1, snf(G)[1].snf) == filter(x -> x != 1, snf(H)[1].snf)
  return b
end

################################################################################
#
#  Direct product
#
################################################################################
#TODO: check the universal properties here!!!
@doc Markdown.doc"""
    direct_product(G::GrpAbFinGen...; task::Symbol = :sum) -> GrpAbFinGen, GrpAbFinGenMap, GrpAbFinGenMap

Returns the direct product $D$ of the abelian groups $G_i$. `task` can be
":sum", ":prod", ":both" or ":none" and determines which canonical maps
are computed as well: ":sum" for the injections, ":prod" for the
projections.
"""
function direct_product(G::GrpAbFinGen...
             ; add_to_lattice::Bool = false, L::GrpAbLattice = GroupLattice, task::Symbol = :sum)
  @assert task in [:prod, :sum, :both, :none]

  Dp = abelian_group(cat([rels(x) for x = G]..., dims = (1,2)))

  set_special(Dp, :direct_product =>G, :show => show_direct_product)
  inj = GrpAbFinGenMap[]
  pro = GrpAbFinGenMap[]
  j = 0
  for g = G
    if task in [:sum, :both]
      m = hom(g, Dp, GrpAbFinGenElem[Dp[j+i] for i = 1:ngens(g)])
      append!(L, m)
      push!(inj, m)
    end
    if task in [:prod, :both]
      m = hom(Dp, g, vcat(GrpAbFinGenElem[g[0] for i = 1:j], gens(g), GrpAbFinGenElem[g[0] for i=j+ngens(g)+1:ngens(Dp)]))
      append!(L, m)
      push!(pro, m)
    end
    j += ngens(g)
  end
  if task == :none
    return Dp
  elseif task == :prod
    return Dp, pro<commit_msg>improve direct_prod by not using hnf<commit_after>################################################################################
#
#       GrpAb/GrpAbFinGen.jl : Finitely generated abelian groups
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: <NAME>, <NAME>
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 <NAME>, <NAME>
#
################################################################################

export abelian_group, free_abelian_group, issnf, ngens, nrels, rels, snf, isfinite,
       isinfinite, rank, order, exponent, istrivial, isisomorphic,
       direct_product, istorsion, torsion_subgroup, sub, quo, iscyclic,
       psylow_subgroup, issubgroup, abelian_groups, flat, tensor_product,
       dual, chain_complex, isexact, homology, free_resolution, obj, map,
       primary_part, isfree

import Base.+, Nemo.snf, Nemo.parent, Base.rand, Nemo.issnf

################################################################################
#
#  Functions for some abstract interfaces
#
################################################################################

elem_type(G::GrpAbFinGen) = GrpAbFinGenElem

elem_type(::Type{GrpAbFinGen}) = GrpAbFinGenElem

parent_type(::Type{GrpAbFinGenElem}) = GrpAbFinGen

isabelian(::GrpAbFinGen) = true

##############################################################################
#
#  Constructors
#
##############################################################################

@doc Markdown.doc"""
    abelian_group(::Type{T} = GrpAbFinGen, M::fmpz_mat) -> GrpAbFinGen

Creates the abelian group with relation matrix `M`. That is, the group will
have `ncols(M)` generators and each row of `M` describes one relation.
"""
abelian_group(M::fmpz_mat; name::String = "") = abelian_group(GrpAbFinGen, M, name=name)

function abelian_group(::Type{GrpAbFinGen}, M::fmpz_mat; name::String = "")
   if issnf(M) && nrows(M) > 0  && ncols(M) > 0 && !isone(M[1, 1])
    N = fmpz[M[i, i] for i = 1:min(nrows(M), ncols(M))]
    if ncols(M) > nrows(M)
      N = vcat(N, fmpz[0 for i = 1:ncols(M)-nrows(M)])
    end
    G = GrpAbFinGen(N)
  else
    G = GrpAbFinGen(M)
  end
  name == "" || set_name!(G, name)
  return G
end

@doc Markdown.doc"""
    abelian_group(::Type{T} = GrpAbFinGen, M::AbstractMatrix{<:IntegerUnion})

Creates the abelian group with relation matrix `M`. That is, the group will
have `ncols(M)` generators and each row of `M` describes one relation.
"""
function abelian_group(M::AbstractMatrix{<:IntegerUnion}; name::String = "")
  return abelian_group(GrpAbFinGen, M, name=name)
end

function abelian_group(::Type{GrpAbFinGen}, M::AbstractMatrix{<:IntegerUnion}; name::String = "")
  return abelian_group(matrix(FlintZZ, M), name=name)
end

function _issnf(N::Vector{T}) where T <: IntegerUnion
  if isone(length(N)) && isone(N[1])
    return false
  end
  for i = 1:length(N)-1
    if isone(abs(N[i]))
      return false
    end
    if iszero(N[i])
      if !iszero(N[i+1])
        return false
      end
    elseif !iszero(mod(N[i+1], N[i]))
      return false
    end
  end
  return true
end

@doc Markdown.doc"""
    abelian_group(::Type{T} = GrpAbFinGen, M::AbstractVector{<:IntegerUnion}) -> GrpAbFinGen
    abelian_group(::Type{T} = GrpAbFinGen, M::IntegerUnion...) -> GrpAbFinGen

Creates the direct product of the cyclic groups $\mathbf{Z}/m_i$,
where $m_i$ is the $i$th entry of `M`.
"""
function abelian_group(M::AbstractVector{<:IntegerUnion}; name::String = "")
  return abelian_group(GrpAbFinGen, M, name=name)
end

function abelian_group(::Type{GrpAbFinGen}, M::AbstractVector{<:IntegerUnion}; name::String = "")
  if _issnf(M)
    G = GrpAbFinGen(M)
  else
	N = zero_matrix(FlintZZ, length(M), length(M))
    for i = 1:length(M)
      N[i,i] = M[i]
    end
    G = GrpAbFinGen(N)
  end
  if !isempty(M)
    res = lcm(M)
    if !iszero(res)
      G.exponent = res
    end
  end
  name == "" || set_name!(G, name)
  return G
end

function abelian_group(M::IntegerUnion...; name::String = "")
  return abelian_group(collect(M), name=name)
end

@doc Markdown.doc"""
    free_abelian_group(::Type{T} = GrpAbFinGen, n::Int) -> GrpAbFinGen

Creates the free abelian group of rank `n`.
"""
free_abelian_group(n::Int) = free_abelian_group(GrpAbFinGen, n)

function free_abelian_group(::Type{GrpAbFinGen}, n::Int)
  return abelian_group(GrpAbFinGen, zeros(Int, n))
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, A::GrpAbFinGen)
  @show_name(io, A)
  @show_special(io, A)

  if issnf(A)
    show_snf(io, A)
  else
    show_gen(io, A)
  end
end

function show_hom(io::IO, G)#::GrpAbFinGen)
  D = get_special(G, :hom)
  D === nothing && error("only for hom")
  print(io, "hom of ")
  print(IOContext(io, :compact => true), D)
end

function show_direct_product(io::IO, G)#::GrpAbFinGen)
  D = get_special(G, :direct_product)
  D === nothing && error("only for direct products")
  print(io, "direct product of ")
  show(IOContext(io, :compact => true), D)
end

function show_tensor_product(io::IO, G)#::GrpAbFinGen)
  D = get_special(G, :tensor_product)
  D === nothing && error("only for tensor products")
  print(io, "tensor product of ")
  show(IOContext(io, :compact => true), D)
end

function show_gen(io::IO, A::GrpAbFinGen)
  print(io, "(General) abelian group with relation matrix\n$(A.rels)")
  if isdefined(A, :snf_map)
    println(io, "\nwith structure of ", domain(A.snf_map))
  end
end

function show_snf(io::IO, A::GrpAbFinGen)
  if get(io, :compact, false)
    print(io, "Abelian group with structure: ")
  else
    print(io, "GrpAb: ")
  end
  show_snf_structure(io, A)
end

function show_snf_structure(io::IO, A::GrpAbFinGen, mul = "x")
  @assert issnf(A)
  len = length(A.snf)
  if len == 0
    print(io, "Z/1")
    return
  end

  if A.snf[1] == 0
    if len == 1
      print(io, "Z")
    else
      print(io, "Z^$(len)")
    end
    return
  end

  i = 1
  while i <= len
    inv = A.snf[i]
    j = 1
    while i + j <= len && inv == A.snf[i + j]
      j += 1
    end
    if iszero(inv)
      print(io, "Z")
      if j > 1
        print(io, "^($j)")
      end
    else
      if j > 1
        print(io, "(Z/$(inv))^$(j)")
      else
        print(io, "Z/$(inv)")
      end
    end
    if i + j - 1 < len
      print(io, " $mul ")
    end
    i += j
  end
end

################################################################################
#
#  Hash function
#
################################################################################

hash_snf(A::GrpAbFinGen, h::UInt) = hash(A.snf, h)

function hash(A::GrpAbFinGen, h::UInt)
  if issnf(A)
    return hash_snf(A, h)
  else
    return hash_snf(snf(A)[1], h)
  end
end

################################################################################
#
#  Field access
#
################################################################################

@doc Markdown.doc"""
    issnf(G::GrpAbFinGen) -> Bool

Returns whether the current relation matrix of the group $G$ is in Smith
normal form.
"""
issnf(A::GrpAbFinGen) = A.issnf

@doc Markdown.doc"""
    ngens(G::GrpAbFinGen) -> Int

Returns the number of generators of $G$ in the current representation.
"""
function ngens(A::GrpAbFinGen)
  if issnf(A)
    return length(A.snf)
  else
    return ncols(A.rels)
  end
end

@doc Markdown.doc"""
    nrels(G::GrpAbFinGen) -> Int

Returns the number of relations of $G$ in the current representation.
"""
function nrels(A::GrpAbFinGen)
  if issnf(A)
    return length(A.snf)
  else
    return nrows(A.rels)
  end
end

@doc Markdown.doc"""
    rels(A::GrpAbFinGen) -> fmpz_mat

Returns the currently used relations of $G$ as a single matrix.
"""
function rels(A::GrpAbFinGen)
  if issnf(A)
    return rels_snf(A)
  else
    return rels_gen(A)
  end
end

function rels_gen(A::GrpAbFinGen)
  return A.rels
end

function rels_snf(A::GrpAbFinGen)
  M = zero_matrix(FlintZZ, ngens(A), ngens(A))
  for i = 1:ngens(A)
    M[i,i] = A.snf[i]
  end
  return M
end

################################################################################
#
#  Hermite normal form
#
################################################################################

function assure_has_hnf(A::GrpAbFinGen)
  if isdefined(A, :hnf)
    return nothing
  end
  R = rels(A)
  if isdefined(A, :exponent) && nrows(R) >= ncols(R)
    A.hnf = hnf_modular_eldiv(R, A.exponent)
  else
    A.hnf = hnf(R)
  end
  return nothing
end

################################################################################
#
#  Smith normal form
#
################################################################################

@doc Markdown.doc"""
    snf(A::GrpAbFinGen) -> GrpAbFinGen, Map

Returns a pair $(G, f)$, where $G$ is an abelian group in canonical Smith
normal form isomorphic to $G$ and an isomorphism $f : G \to A$.
"""
function snf(G::GrpAbFinGen)
  if isdefined(G, :snf_map)
    return domain(G.snf_map)::GrpAbFinGen, G.snf_map::GrpAbFinGenMap
  end

  if issnf(G)
    G.snf_map = GrpAbFinGenMap(G) # identity
    return G, G.snf_map::GrpAbFinGenMap
  end
  if isdefined(G, :exponent)
    if isdefined(G, :hnf)
      S, T = snf_for_groups(G.hnf, G.exponent)
    else
      S, T = snf_for_groups(G.rels, G.exponent)
    end
  else
    S, _, T = snf_with_transform(G.rels, false, true)
  end

  m = min(nrows(S), ncols(S))
  if m > 0 && nrows(S) >= ncols(S)
    e = S[m, m]
    if e > 1
      if fits(Int, e) && isprime(e)
        F = GF(Int(e), cached = false)
        TF = map_entries(F, T)
        iT = lift(inv(TF))
      else
        iT = invmod(T, e)
      end
    else
      iT = inv(T)
    end
  else
    iT = inv(T)
  end
  return _reduce_snf(G, S, T, iT)
end

# For S in SNF with G.rels = U*S*T and Ti = inv(T) this removes
# the ones at the diagonal of S and constructs the homomorphism.
function _reduce_snf(G::GrpAbFinGen, S::fmpz_mat, T::fmpz_mat, Ti::fmpz_mat)

  d = fmpz[S[i,i] for i = 1:min(nrows(S), ncols(S))]

  while length(d) < ngens(G)
    push!(d, 0)
  end

  pos = Int[i for i = 1:length(d) if !isone(d[i])]
  r = Int[i for i = 1:nrows(T)]
  s = fmpz[ d[i] for i in pos]
  TT = sub(T, r, pos)
  TTi = sub(Ti, pos, r)

  H = GrpAbFinGen(s)
  if !isempty(s) && !iszero(s[end])
    H.exponent = s[end]
    G.exponent = s[end]
  end
  mp = hom(H, G, TTi, TT, check = false)
  G.snf_map = mp
  return H, mp::GrpAbFinGenMap
end

################################################################################
#
#  Finiteness
#
################################################################################

@doc Markdown.doc"""
    isfinite(A::GrpAbFinGen) -> Bool

Returns whether $A$ is finite.
"""
isfinite(A::GrpAbFinGen) = issnf(A) ? isfinite_snf(A) : isfinite_gen(A)

isfinite_snf(A::GrpAbFinGen) = length(A.snf) == 0 || !iszero(A.snf[end])

isfinite_gen(A::GrpAbFinGen) = isfinite(snf(A)[1])

@doc Markdown.doc"""
    isinfinite(A::GrpAbFinGen) -> Bool

Returns whether $A$ is infinite.
"""
isinfinite(A::GrpAbFinGen) = !isfinite(A)

################################################################################
#
#  Rank
#
################################################################################

@doc Markdown.doc"""
    rank(A::GrpAbFinGen) -> Int

Returns the rank of $A$, that is, the dimension of the
$\mathbf{Q}$-vectorspace $A \otimes_{\mathbf Z} \mathbf Q$.
"""
rank(A::GrpAbFinGen) = issnf(A) ? rank_snf(A) : rank_gen(A)

rank_snf(A::GrpAbFinGen) = length(findall(x -> iszero(x), A.snf))

rank_gen(A::GrpAbFinGen) = rank(snf(A)[1])

rank(A::GrpAbFinGen, p::Union{Int, fmpz}) = issnf(A) ? rank_snf(A, p) : rank_snf(snf(A)[1], p)

function rank_snf(A::GrpAbFinGen, p::Union{Int, fmpz})
  if isempty(A.snf)
    return 0
  end
  if !iszero(mod(A.snf[end], p))
    return 0
  end
  i = findfirst(x -> iszero(mod(x, p)), A.snf)
  return length(A.snf)-i+1
end


################################################################################
#
#  Order
#
################################################################################

@doc Markdown.doc"""
    order(A::GrpAbFinGen) -> fmpz

Returns the order of $A$. It is assumed that $A$ is finite.
"""
order(A::GrpAbFinGen) = issnf(A) ? order_snf(A) : order_gen(A)

function order_snf(A::GrpAbFinGen)
  isinfinite(A) && error("Group must be finite")
  return prod(A.snf)
end

order_gen(A::GrpAbFinGen) = order(snf(A)[1])

################################################################################
#
#  Exponent
#
################################################################################

@doc Markdown.doc"""
    exponent(A::GrpAbFinGen) -> fmpz

Returns the exponent of $A$. It is assumed that $A$ is finite.
"""
function exponent(A::GrpAbFinGen)
  if issnf(A)
    res = exponent_snf(A)
    if !iszero(res)
      A.exponent = res
    end
    return res
  else
    res = exponent_gen(A)
    if !iszero(res)
      A.exponent = res
    end
    return res
  end
end

function exponent_snf(A::GrpAbFinGen)
  isinfinite(A) && error("Group must be finite")
  ngens(A)==0 && return fmpz(1)
  return A.snf[end]
end

exponent_gen(A::GrpAbFinGen) = exponent(snf(A)[1])

################################################################################
#
#  Trivial
#
################################################################################

@doc Markdown.doc"""
    istrivial(A::GrpAbFinGen) -> Bool

Checks if $A$ is the trivial group.
"""
istrivial(A::GrpAbFinGen) = isfinite(A) && isone(order(A))

################################################################################
#
#  Isomorphism
#
################################################################################

@doc Markdown.doc"""
    isisomorphic(G::GrpAbFinGen, H::GrpAbFinGen) -> Bool

Checks if $G$ and $H$ are isomorphic.
"""
function isisomorphic(G::GrpAbFinGen, H::GrpAbFinGen)
  b = filter(x -> x != 1, snf(G)[1].snf) == filter(x -> x != 1, snf(H)[1].snf)
  return b
end

################################################################################
#
#  Direct product
#
################################################################################
#TODO: check the universal properties here!!!
@doc Markdown.doc"""
    direct_product(G::GrpAbFinGen...; task::Symbol = :sum) -> GrpAbFinGen, GrpAbFinGenMap, GrpAbFinGenMap

Returns the direct product $D$ of the abelian groups $G_i$. `task` can be
":sum", ":prod", ":both" or ":none" and determines which canonical maps
are computed as well: ":sum" for the injections, ":prod" for the
projections.
"""
function direct_product(G::GrpAbFinGen...
             ; add_to_lattice::Bool = false, L::GrpAbLattice = GroupLattice, task::Symbol = :sum)
  @assert task in [:prod, :sum, :both, :none]

  Dp = abelian_group(cat([rels(x) for x = G]..., dims = (1,2)))
  for x = G
    assure_has_hnf(x)
  end
  Dp.hnf = cat([x.hnf for x = G]..., dims = (1,2))

  set_special(Dp, :direct_product =>G, :show => show_direct_product)
  inj = GrpAbFinGenMap[]
  pro = GrpAbFinGenMap[]
  j = 0
  for g = G
    if task in [:sum, :both]
      m = hom(g, Dp, GrpAbFinGenElem[Dp[j+i] for i = 1:ngens(g)], check = false)
      append!(L, m)
      push!(inj, m)
    end
    if task in [:prod, :both]
      m = hom(Dp, g, vcat(GrpAbFinGenElem[g[0] for i = 1:j], gens(g), GrpAbFinGenElem[g[0] for i=j+ngens(g)+1:ngens(Dp)]), check = false)
      append!(L, m)
      push!(pro, m)
    end
    j += ngens(g)
  end
  if task == :none
    return Dp
  elseif task == :prod
    return Dp, pro
{'fraction_non_alphanumeric': 0.2689985052316891, 'fraction_numerical': 0.004783258594917788, 'mean_word_length': 4.430357142857143, 'pattern_counts': {'":': 12, '<': 33, '<?xml version=': 0, '>': 81, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5685860', 'n_tokens_mistral': 12084, 'n_tokens_neox': 11521, 'n_words': 2970}