################################################################################
#
#  Two element generated ideals
#
################################################################################

@doc Markdown.doc"""
    has_2_elem(A::NfAbsOrdIdl) -> Bool

Returns whether $A$ is generated by two elements.
"""
function has_2_elem(A::NfAbsOrdIdl)
  return isdefined(A, :gen_two)
end

@doc Markdown.doc"""
    has_weakly_normal(A::NfAbsOrdIdl) -> Bool

Returns whether $A$ has weakly normal two element generators.
"""
function has_weakly_normal(A::NfAbsOrdIdl)
  return (isdefined(A, :gens_weakly_normal) &&
        A.gens_weakly_normal == true) || has_2_elem_normal(A)
end

@doc Markdown.doc"""
    has_2_elem_normal(A::NfAbsOrdIdl) -> Bool

Returns whether $A$ has normal two element generators.
"""
function has_2_elem_normal(A::NfAbsOrdIdl)
  #the one ideal <1, ?> is automatomatically normal>
  return isdefined(A, :gens_normal) && (A.gen_one == 1 || A.gens_normal > 1)
end

################################################################################
#
#  Predicates
#
################################################################################

function isconsistent(A::NfAbsOrdIdl)
  if has_2_elem_normal(A) && !iszero(A.gen_two)
    if !defines_2_normal(A) 
      return false
    end
  end
  if has_norm(A)
    b = basis_mat(A, copy = false)
    if det(b) != A.norm
      return false
    end
  end
  return true
end

# check if gen_one,gen_two is a P(gen_one)-normal presentation
# see Pohst-Zassenhaus p. 404
function defines_2_normal(A::NfAbsOrdIdl)
  m = A.gen_one
  gen = A.gen_two
  mg = denominator(inv(gen.elem_in_nf), order(A))
  # the minimum of ideal generated by g
  g = gcd(m,mg)
  return gcd(m, div(mg, g)) == 1
end

###########################################################################################
#
#  2-element normal presentation
#
###########################################################################################

# The following makes sure that A has a weakly normal presentation
# Recall that (x,y) are a weakly normal presentation for A
# if and only if norm(A) = gcd(norm(x), norm(y))
#
# Maybe we should allow an optional paramter (an fmpz),
# which should be the first generator.
# So far, the algorithm just samples (lifts of) random elements of A/m^2,
# where m is the minimum of A.

function _assure_weakly_normal_presentation(A::NfAbsOrdIdl)
  if has_2_elem(A) && has_weakly_normal(A)
    return nothing
  end

  if isdefined(A, :princ_gen)
    x = A.princ_gen
    b = x.elem_in_nf

    bi = inv(b)

    A.gen_one = denominator(bi, order(A))
    A.minimum = A.gen_one
    A.gen_two = x
    if !isdefined(A, :norm)
      A.norm = abs(_normmod(A.gen_one^degree(order(A)), b))
    end
    @hassert :NfOrd 1 gcd(A.gen_one^degree(order(A)),
                    FlintZZ(norm(A.gen_two))) == A.norm

    if A.gen_one == 1
      A.gens_normal = 2*A.gen_one
    else
      A.gens_normal = A.gen_one
    end
    A.gens_weakly_normal = 1
    return nothing
  end

  O = order(A)
  if minimum(A) == 0
    A.gen_one = minimum(A)
    A.gen_two = zero(O)
    A.gens_weakly_normal = true
    return nothing
  end

  if minimum(A) == 1
    A.gen_one = minimum(A)
    A.gen_two = one(O)
    A.gens_weakly_normal = true
    A.gens_normal = fmpz(2)
    return nothing
  end

  if !issimple(nf(order(A))) || !isdefining_polynomial_nice(nf(order(A)))
    B = Array{fmpz}(undef, degree(O))
    Amin2 = minimum(A)^2
    Amind = minimum(A)^degree(O)
    BB = basis(A, copy = false)
    r = -Amin2:Amin2

    while true
      rand!(B, r)
      gen = dot(B, BB)
      if norm(A) == gcd(Amind, norm(gen))
        A.gen_one = minimum(A)
        A.gen_two = O(gen, false)
        A.gens_weakly_normal = true
        return nothing
      end
    end
  end

  Amin2 = minimum(A)^2
  Amind = gcd(minimum(A)^degree(O), minimum(A)*norm(A))

  B = Array{fmpz}(undef, degree(O))

  gen = O()

  r = -Amin2:Amin2

  m = zero_matrix(FlintZZ, 1, degree(O))

  cnt = 0
  while true
    cnt += 1

    if cnt > 100 && is2_normal_difficult(A)
      assure_2_normal_difficult(A)
      return nothing
    end

    #if cnt > 1000
    #  println("Having a hard time find weak generators for $A")
    #end

    rand!(B, r)

    # Put the entries of B into the (1 x d)-Matrix m
    for i in 1:degree(O)
      s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), m, 0, i - 1)
      ccall((:fmpz_set, :libflint), Nothing, (Ptr{fmpz}, Ref{fmpz}), s, B[i])
    end

    if iszero(m)
      continue
    end

    mul!(m, m, basis_mat(A, copy = false))
    d = denominator(basis_mat(O, copy = false))
    mul!(m, m, basis_mat(O, copy = false).num)
    gen = elem_from_mat_row(nf(O), m, 1, d)
    d = denominator(gen)
    f, e = ppio(d, minimum(A, copy = false))
    gen = mod(numerator(gen), f*minimum(A)^2)//f
    if iszero(gen)
      continue
    end

    # the following should be done inplace
    #gen = dot(reshape(Array(mm), degree(O)), basis(O))
    if norm(A) == _normmod(Amind, O(gen, false))#gcd(Amind, numerator(norm(gen)))
      A.gen_one = minimum(A)
      A.gen_two = O(gen, false)
      A.gens_weakly_normal = true
      return nothing
    end
  end
end
#
# Here some random guesses for the difficult 2-element generators
# degree  | d
#   < 7   | 1
#  8 - 12 | 2 * 3
# 13 - 20 | 2 * 3 * 5
#  >= 21  | 2 * 3 * 5 * 7

function is2_normal_difficult(A::NfAbsOrdIdl)
  d = 2
  m = minimum(A)
  ZK = order(A)
  n = degree(ZK)

  if n < 7
    return false
  end

  if n < 12 && isone(gcd(m, 2 * 3))
    return false
  end

  if n < 20 && isone(gcd(m, 2 * 3 * 5))
    return false
  end

  if isone(gcd(m, 2 * 3 * 5 * 7))
    return false
  end

  return true
end

function valuation_of_ideal_difficult(A::NfAbsOrdIdl, P::NfAbsOrdIdl)
  B = basis(A, copy = false)
  res = valuation(B[1], P)
  for i = 2:length(B)
    if iszero(res)
      break
    end
    res = min(res, valuation(B[i], P))
  end
  return res
end


function assure_2_normal_difficult(A::NfAbsOrdIdl)

  m = minimum(A)
  ZK = order(A)
  n = degree(ZK)

  if !is2_normal_difficult(A)
    assure_2_normal(A)
    return nothing
  end

  if n < 12
    d = Int[2, 3]
  elseif n < 20
    d = Int[2, 3, 5]
  else
    d = Int[2, 3, 5, 7]
  end
  
  I = ideal(ZK, 1)
  for i = 1:length(d)
    m1, m = ppio(m, fmpz(d[i]))
    if isone(m1)
      continue
    end
    A1 = gcd(A, m1)
    lp = prime_decomposition(ZK, d[i])
    v = Int[valuation_of_ideal_difficult(A1, p[1]) for p = lp]
    for i = 1:length(v)
      if v[i] > 0 
        I *= lp[i][1]^v[i]
      end
    end
  end
  #Now the ideal I has a 2 normal presentation
  if isone(m)
    A.gen_one = I.gen_one
    A.gen_two = I.gen_two
    A.gens_normal = I.gens_normal
    A.gens_weakly_normal = I.gens_weakly_normal
    A.gens_short = I.gens_short
    return nothing
  end
  J = gcd(A, m)
  assure_2_normal(J)
  C = I * J
  A.gen_one = C.gen_one
  A.gen_two = C.gen_two
  A.gens_normal = C.gens_normal
  A.gens_weakly_normal = C.gens_weakly_normal
  A.gens_short = C.gens_short
  return nothing
end

################################################################################
#
#  Assure 2 normal
#
################################################################################

function assure_2_normal(A::NfAbsOrdIdl)
  if has_2_elem(A) && has_2_elem_normal(A)
    @hassert :NfOrd 1 isconsistent(A)
    return nothing
  end
  O = order(A)
  K = nf(O)
  n = degree(K)

  if norm(A) == 1
    A.gen_one = fmpz(1)
    A.gen_two = one(O)
    A.gens_normal = fmpz(1)
    return
  end

  if has_2_elem(A)
    m = minimum(A)
    bas = basis(O)
    # Magic constants
    if m > 1000
      r = -500:500
    else
      r = -div(Int(m)+1,2):div(Int(m)+1,2)
    end
    #gen = K()
    #s = K()
    gen = zero(O)
    s = O()
    cnt = 0
    while true
      cnt += 1
      if cnt > 100 && is2_normal_difficult(A)
        assure_2_normal_difficult(A)
        @hassert :NfOrd 1 isconsistent(A)
        return  
      end
      if cnt > 1000
        error("Having a hard time making generators normal for $A")
      end
      rand!(s, O, r)
      mul!(s, s, A.gen_two)
      add!(gen, rand(r)*A.gen_one, gen)
      add!(gen, s, gen)
      gen = mod(gen, m^2)

      if iszero(gen)
        continue
      end

      mg = _minmod(m^2, gen)
      g = gcd(m, mg)
      if gcd(m, div(mg, g)) == 1
        if gcd(m^n, norm(gen)) != norm(A)
          @vprint :NfOrd 2 "\n\noffending ideal $A \ngen is $gen\nWrong ideal\n"
          cnt += 10
          continue
        end
        break
      end
    end
    @vprint :NfOrd 2 "used $cnt attempts\n"
    A.gen_one = m
    A.gen_two = gen
    A.gens_normal = m
    @hassert :NfOrd 1 isconsistent(A)
    return
  end

  m = minimum(A)
  bas = nf_elem[ elem_in_nf(b) for b in basis(A, copy = false)]
  # Magic constants
  if m > 1000
    r = -500:500
  else
    r = -div(Int(m)+1,2):div(Int(m)+1,2)
  end
  cnt = 0
  while true
    cnt += 1
    if cnt > 100 && is2_normal_difficult(A)
      assure_2_normal_difficult(A)
      @hassert :NfOrd 1 isconsistent(A)
      return  
    end
    if cnt > 1000
      error("Having a hard time making generators normal for $A")
    end
    gen = O(rand(bas, r))
    gen = mod(gen, m^2)

    if iszero(gen)
      continue
    end

    mg = _minmod(m^2, gen)
    if gcd(m, div(mg, gcd(mg, m))) == 1
      if gcd(m^n, norm(gen)) != norm(A)
        @vprint :NfOrd 1 "\n\noffending ideal $A \ngen is $gen\nWrong ideal\n"
        cnt += 10
        continue
      end
      break
    end
  end
  @vprint :NfOrd 2 "used $cnt attempts\n"
  A.gen_one = m
  A.gen_two = gen
  A.gens_normal = m
  @hassert :NfOrd 1 isconsistent(A)
  return nothing
end

