text:
# Genus representatives for quadratic lattices,
#
# With permission ported from Magma package of Markus Kirschmer:
# http://www.math.rwth-aachen.de/~Markus.Kirschmer/magma/lat.html

add_verbose_scope(:GenRep)
add_assert_scope(:GenRep)

export genus_representatives

################################################################################
#
#  SpinorGeneraCtx
#
################################################################################

# To keep track of ray class groups
mutable struct SpinorGeneraCtx
  mR::MapRayClassGrp # ray class group map
  mQ::GrpAbFinGenMap # quotient
  rayprimes::Vector{NfAbsOrdIdl{AnticNumberField, nf_elem}}
  criticalprimes::Vector{NfAbsOrdIdl{AnticNumberField, nf_elem}}

  function SpinorGeneraCtx()
    return new()
  end
end

function SpinorGeneraCtx(L::QuadLat)
  R = base_ring(L)
  F = number_field(R)

  RCG, mRCG, Gens = _compute_ray_class_group(L)

  # 1) Map the generators into the class group to create the factor group.

  subgroupgens = GrpAbFinGenElem[_map_idele_into_class_group(mRCG, [g]) for g in Gens ]

  for g in gens(RCG)
    push!(subgroupgens, 2*g)
  end

  Q, mQ = quo(domain(mRCG), subgroupgens)

  @vprint :GenRep 1 "Ray class group: size = $(order(RCG))\n"
  @vprint :GenRep 1 "Ray class group quotient: size = $(order(Q))
  (this is the number of spinor + genera in Genus(L))\n"

  # 2) Find good generators (that is: generators which are not dyadic, and not
  #    in BadPrimes(L) -- so that neighbour generation is always possible),
  #    of this factor group.

  #   Only pick ideles of the form (1,...,1,p,1,...,1) so that:
  #   a) nothing has to be corrected before they can be mapped down into the factor group
  #   b) we can simply store the prime ideal and know that it corresponds to
  #      the (1,...,1,p,1,...,1) These prime ideals will be stored in
  #      SpinorGeneraCtx.criticalprimes.
  #
  #   the primes which generate the spinor genera in Genus(L):
  #   if L is not definite, we need one prime for each factor group element.
  #   if L is definite, we only need a generating set of primes for the factor group.

  inf_plc = defining_modulus(mRCG)[2]

  critical_primes = _get_critical_primes(L, mRCG, inf_plc, mQ, true) # fulllist = not isdefinite?
  @vprint :GenRep 1 "good primes over $([minimum(q) for q in critical_primes])
  (together with the squares) generate the subgroup.\n"

  res = SpinorGeneraCtx()

  res.mR = mRCG
  res.mQ = mQ
  res.criticalprimes = critical_primes
  res.rayprimes = collect(keys(mRCG.fact_mod))

  return res
end

################################################################################
#
#  Genus representatives for quadratic lattices
#
################################################################################

@doc Markdown.doc"""
    genus_representatives(L::QuadLat; max = inf, use_auto = false)
                                                        -> Vector{QuadLat}

Computes representatives for the isometry classes in the genus of $L$.

At most `max` representatives are returned. The use of automorphims can be
disabled by
"""
function genus_representatives(L::QuadLat; max = inf, use_auto = true, use_mass = true)
  # Otherwise the isomorphism to the class group fails, cf. §102 in O'Meara.
  @req max >= 1 "Must find at least one representative"

  if rank(L) == 1
    return typeof(L)[L]
  end

  if rank(L) == 2
    if isdefinite(L)
      return _genus_representatives_binary_quadratic_definite(L, max = max, use_auto = true, use_mass = true)
    else
      @req degree(base_ring(L)) == 1 "Binary indefinite quadratic lattices must be only over the rationals"
      return _genus_representatives_binary_quadratic_indefinite(L)
    end
  end

  if !isdefinite(L)
    @vprint :GenRep 1 "Genus representatives of indefinite lattice\n"
    return spinor_genera_in_genus(L, [])
  end

  res = typeof(L)[]

  p = _smallest_norm_good_prime(L)

  spinor_genera = spinor_genera_in_genus(L, typeof(p)[p])

  if use_mass
    @vprint :GenRep 1 "Computing mass exactly ...\n"
    _mass = mass(L)
    @vprint :GenRep 1 "... $(_mass)\n"
  else
    _mass = -one(FlintQQ)
  end

  @vprint :GenRep 1 "Found $(length(spinor_genera)) many spinor genera in genus\n"

  for LL in spinor_genera
    @hassert :GenRep 3 all(!isisometric(X, LL)[1] for X in res)
    new_lat =  iterated_neighbours(LL, p, use_auto = use_auto,
                                          max = max - length(res),
                                          mass = _mass//length(spinor_genera))
    append!(res, new_lat)
  end

  if max > length(res) && use_mass
    if sum(fmpq[1//automorphism_group_order(LL) for LL in res]) != _mass
      throw(error("Something very wrong"))
    end
  end

  return res
end

function spinor_genera_in_genus(L, mod_out)
#{A sequence of lattices representing the spinor genera in the genus of L}
  @req rank(L) >= 3 "(Currently) rank must be at least 3"
  # Otherwise the isomorphism to the class group fails, cf. §102 in O'Meara.

  # We have to find out whether ProperSpinorGenus == SpinorGenus.

  # 1) Find a nonzero element in the norm group of L.
  C = SpinorGeneraCtx(L)

  # We don't need minimal generators, to make them not too large.
  Gr = gram_matrix_of_generators(L, minimal = false)

  R = base_ring(L)
  F = nf(R)

  spinornorm = zero(F)

  # The smaller the element, the better
  for d in diagonal(Gr)
    if iszero(spinornorm) || (abs(norm(d)) < abs(norm(spinornorm)))
      spinornorm = d
    end
  end
  if iszero(spinornorm)
    i = 1
    while iszero(Gr[1, i])
      i += 1
      if i > ncols(Gr)
        throw(error("Lattice is degenerated"))
      end
      @assert !iszero(Gr[1, i])
      spinornorm = 2 * Gr[1, i]
    end
  end

  # 2) At a place p where spinornorm does not generate norm(L_p)
  #    we add <p, spinornorm * normgenerator > to the idele

  nor = norm(L)

  differ = ideal_type(R)[]
  for p in vcat(support(nor), support(spinornorm * R))
    if valuation(nor, p) != valuation(spinornorm, p)
      push!(differ, p)
    end
  end

  idele = Tuple{ideal_type(R), elem_type(F)}[]

  for p in unique(vcat(differ, C.rayprimes))
    local _norm::elem_type(F)
    if isdyadic(p)
      _,_,_,a = _genus_symbol_kirschmer(L, p).data
      _norm = a[1]::elem_type(F)
    else
      GS = _genus_symbol_kirschmer(L, p)
      G = GS.data
      pi = GS.x::elem_type(F)
      _norm = pi^(G[1][2]::Int)
      if G[1][1]::Int == 1 && G[1][3]::Int == -1
        k, h = ResidueField(R, p)
        _norm = _norm * elem_in_nf(h\non_square(k))
      end
    end
    push!(idele, (p, spinornorm * _norm))
  end

  # 3) Map the idele into the idele class group.
  #    Then h == 0 iff SpinorGenus == ProperSpinorGenus

  g = _map_idele_into_class_group(C.mR, idele)
  h = C.mQ(g)

  # 4) Now there are 2^#primes spinor genera, which can be reached by neighbours:

  modout = push!(elem_type(codomain(C.mQ))[ C.mQ(_map_idele_into_class_group(C.mR,  [(p, elem_in_nf(uniformizer(p)))])) for p in mod_out], h)
  primes = _spinor_generators(L, C, modout)

  res = typeof(L)[ L ]

  for p in primes
    N = neighbours(L, p, max = 1)[1]
    pN = pseudo_matrix(N)
    for i in 1:length(res)
      pM = pseudo_matrix(res[i])
      new_pm = _sum_modules(_module_scale_ideal(p, pM), pN)
      new_pm = _intersect_modules(new_pm, _module_scale_ideal(inv(p), pM))
      LL = lattice(ambient_space(L), new_pm)
      # LL = (p*M + N) \cap p^-1 * M
      push!(res, LL)
    end
  end

  return res
end

function _smallest_norm_good_prime(L)
  OK = base_ring(L)
  lp = ideal_type(OK)[p for p in bad_primes(L, even = true) if isdyadic(p) || !ismodular(L, p)[1]]
  limit = 20
  while true
    lq = prime_ideals_up_to(OK, limit)
    for q in lq
      if !(q in lp)
        return q
      end
    end
    limit = 2 * limit
    if limit > 2^8
      throw(error("Something off"))
    end
  end
end

# The spinor norm of L at p, as calculated by C. Beli for dyadic p and by M.
# Kneser for non-dyadic p.
#
# Returns a subspace of local_multiplicative_group_modulo_squares(p), and a
# boolean which is true iff the spinor norm is exactly the units
function spinor_norm(L, p)
  V, g = local_multiplicative_group_modulo_squares(p)
  R = base_ring(L)
  e = valuation(R(2), p)
  if !isdyadic(p)
    # cf. Satz 3 in Kneser, "Klassenzahlen indefiniter quadratischer Formen in
    # drei oder mehr Veränderlichen":
    J, G, E = jordan_decomposition(L, p)
    if any(g -> ncols(g) >= 2, G)
      @vprint(:GenRep, 1,"""Spinor norm over $(minimum(p))
  This lattice has a 2-dimensional Jordan constituent, and p is odd. Spinor norm
  is either F^* or O_F^*(F^*)^2, i.e. we will find a vector space of dimension
  $(gens(V)) or $(ngens(V) - 1).\n""")
      # Which of the two is the case?
      # TODO: It is not a good idea to rely on implementation details of
      #       local_multiplicative_group_modulo_squares
      if length(unique([e % 2 for e in E])) == 1
        return gens(V)[1:ngens(V) - 1], V, g, true
      else
        return gens(V), V, g, false
      end
    end
    # This is the obscure case where p is odd and each of its Jordan components
    # has dimension 1. In this case, the units need not be contained in the
    # Spinor norm.
    #
    # Generators of the (principal) norm ideals of the Jordan components: since
    # p is odd, the norms (and the scales) are just the dimensions of the
    # Jordan components
     normgens = nf_elem[g[1,1] for g in G]

     twonormgens = nf_elem[]

     for i in 1:length(normgens)
       for j in 1:length(normgens)
         push!(twonormgens, normgens[i] * normgens[j])
       end
     end

     twonormvectors = [g\(x) for x in twonormgens]

     @vprint :GenRep "Spinor norm odd p, norm generators of the $(length(G)) Jordan components are: $(normgens), $(twonormgens) $(twonormvectors)"
    # cf. Kneser 1956, Satz 3:
    _SN, mS = sub(V, twonormvectors)
    #@assert length(rels(_SN)) == 0 # free
    SN = [ mS(s) for s in gens(_SN) ]
  else
    bong = good_bong(L, p)
    @hassert :GenRep 1 isgood_bong(bong, p)
    if !has_propertyA(L, p)
      @vprint(:GenRep, 1,"""Spinor norm over dyadic prime:
  This lattice does not have property A. Spinor norm is either F^* or
  O_F^*(F^*)^2, i.e. we will find a vector space of dimension $(ngens(V)) or
  $(ngens(V) - 1)\n""")
      # Using section 7, Thm. 3 in Beli 2002, one can decide which of the two
      # cases applies. This is why we needed to compute a *good* BONG:
      for i in 1:(length(bong) - 1)
        BG, mBG = G_function(bong[i + 1]//bong[i], V, g, p)
        for bg in gens(BG)
          if mBG(bg)[ngens(V)] != 0
            # the whole group
            return gens(V), V, g, false
          end
        end
      end

      for i in 1:(length(bong) - 2)
        if valuation(bong[i], p) == valuation(bong[i + 2], p)
          @assert iszero(mod(valuation(bong[i + 1], p) - valuation(bong[i], p), 2))
          if mod(div(valuation(bong[i + 1], p) - valuation(bong[i], p), 2), 2) != mod(e, 2)
            # the whole group
            return gens(V), V, g, false
          end
        end
      end

      # If all checks have passed, the spinor norm is exactly the units:
      return gens(V)[1:(ngens(V) - 1)], V, g, true
    end
    # cf. Beli 2003, Thm. 1
    SN = elem_type(V)[]::Vector{elem_type(V)}
    for i in 2:rank(L)
      #SN = SN + G_function(bong[i]//bong[i - 1], V, g, p)
      _G, _mG = G_function(bong[i]//bong[i - 1], V, g, p)
      new_gens = append!(SN, elem_type(V)[_mG(g) for g in gens(_G)])
      _SN, mS = sub(V, new_gens)
      SN = elem_type(V)[ mS(s) for s in gens(_SN) ]
    end
    # For why we can take the Minimum in what follows here, see the remark on p. 161 in Beli 2003:
    k = findfirst(i -> mod(valuation(bong[i + 2], p) - valuation(bong[i], p), 2) == 0, 1:(rank(L) - 2))
    if k !== nothing
      alpha = minimum(div(valuation(bong[i + 2], p) - valuation(bong[i], p), 2) for i in 1:(rank(L) - 2) if mod(valuation(bong[i + 2], p) - valuation(bong[i], p), 2) == 0)
      # SN = SN + one_plus_power_of_P(alpha, V, g, P)
      _G, _mG = _one_plus_power_of_p(alpha, V, g, p)
      _SN, mS = sub(V, append!(SN, [_mG(g) for g in gens(_G)]))
      #@assert length(rels(_SN)) == 0 # free
      SN = elem_type(V)[ mS(s) for s in gens(_SN) ]
    end
  end
  # Test if S is equal to SN
  S, mS = sub(V, gens(V)[1:(ngens(V) - 1)])
  #@assert length(rels(S)) == 0
  W,_ = sub(V, append!(gens(V)[1:(ngens(V) - 1)], SN))
  #@assert length(rels(W)) == 0
  if length(SN) == ngens(S) && ngens(W) == ngens(V) - 1 # this means SN + V[1:dim(V) -1] == V[1:dim(V) - 1]
    fl = true
  else
    fl = false
  end

  return SN, V, g, fl
end

# Return a good BONG of L at a dyadic prime p, as defined by C. Beli.
function good_bong(L, p)
# Any BONG obtained from a maximal norm splitting is good, cf. Corollary 4.2 in Beli.
# If a maximal norm splitting is calculated, a good BONG can be read off by joining
# together the BONGs of the 1- and 2-dimensional components.
  @req base_ring(L) == order(p) "Incompatible arguments"
  @req isdyadic(p) "Prime must be dyadic"
  G, JJ = maximal_norm_splitting(L, p)
  K = nf(base_ring(L))
  bong = nf_elem[]
  for i in 1:length(G)
    GG = G[i]
    if nrows(GG) == 2
      bong = append!(bong, _make_bong_dim_2(quadratic_lattice(K, identity_matrix(K, 2), gram_ambient_space = GG), p))
    elseif nrows(GG) == 1
      push!(bong, GG[1, 1])
    else
      throw(error("This should not happen"))
    end
  end
  return bong
end

# A maximal norm splitting of L at a dyadic prime p, as defined by C. Beli.
# Returns a Jordan decomposition into 1x1 and 2x2 components, and the
# corresponding list of basis vectors.
function maximal_norm_splitting(L, p)
  # Follow Beli, "Representation of Integral Quadratic Forms over Dyadic Local Fields", section 7
  @req base_ring(L) == order(p) "Incompatible arguments"
  @req isdyadic(p) "Prime must be dyadic"

  R = base_ring(L)
  K = nf(R)
  e = ramification_index(p)
  uni = uniformizer(p)
  J, G, _ = jordan_decomposition(L, p)
  # join J into one matrix of base vectors
  JJ = reduce(vcat, J)
  # join the individual Gram matrices:
  A = diagonal_matrix(G)

  # read the finer decomposition:
  G = dense_matrix_type(elem_type(K))[]
  steps = UnitRange{Int}[]
  i = 1
  n = ncols(A)
  while i <= n
    if i != n && !iszero(A[i, i + 1])
      size = 2
    else
      size = 1
    end
    push!(G, sub(A, i:(i + size - 1), i:(i + size - 1)))
    push!(steps, i:(i + size - 1))
    i = i + size
  end

  if all(s -> length(s) == 1, steps)
    return G, JJ
  end

  # Now turn this decomposition into unary/binary components into a maximal norm splitting
  failset = Int[]
  while true
    #  Scale valuations sL, norm generators aL, norm valuations uL, weight
    #  valuations wL of the decomposition
    sL, aL, uL, wL = _scales_and_norms(G, p, uni)

    failset_old = failset
    b, i, failset = __ismaximal_norm_splitting(G, sL, aL, p)
    @assert isempty(failset_old) || length(failset_old) > length(failset)
    if b
      break
    end
    @assert i >= 0
    # If b is false and i == 0, something is wrong with our lattice

    # Here comes the interesting stuff:
    @assert length(steps[i]) == 2
    # Unary components already have maximal norm.

    # The maximal norm splitting condition is violated at index i.
    find_j = append!([valuation(aL[k], p) for k in (i+1):length(aL)], [2*(sL[i] - sL[k]) + valuation(aL[k], p) for k in 1:(i - 1)])
    @assert length(find_j) > 0
    min, j = findmin(find_j)
    if j <= length(aL) - i
      j = j + i
      # We want j to correspond to a Jordan component, not an index in find_j
      @assert j > i
      # This is case 1 in Beli.
      beli_correction!(L, G, JJ, steps, i, j, p)
    else
      j = j - (length(aL) - i)
      # We want j to correspond to a Jordan component, not an index in find_j
      @assert j < i
      # This is case 2 in Beli.
      # Switch to a basis of the dual lattice L^#:
      for k in 1:length(steps)
        for l in 1:length(steps[k])
          for u in 1:ncols(JJ)
            JJ[steps[k][l], u] = JJ[steps[k][l], u] * uni^(-sL[k])
          end
        end
        B = sub(JJ, steps[k][1]:steps[k][length(steps[k])], 1:ncols(JJ))

        @assert valuation(scale(quadratic_lattice(K, identity_matrix(K, nrows(B)), gram_ambient_space = B * gram_matrix(ambient_space(L)) * B')), p) == -sL[k]
      end
      # Apply case 1 to the reversed orthogonal sum of the dual lattices:

      steps = reverse!(steps)
      # Update Gram matrices for the reversed, dualized lattice:
      for k in 1:length(G)
        B = sub(JJ, steps[k][1]:steps[k][length(steps[k])], 1:ncols(JJ))
        G[k] = B * gram_matrix(ambient_space(L)) * B'
        @assert nrows(G[k]) in [1, 2]
      end
      # Component i is now at position #aL-i+1
      @assert length(steps[length(aL) - i + 1]) == 2
      beli_correction!(L, G, JJ, steps, length(aL) - i + 1, length(aL) - j + 1, p)
      steps = reverse!(steps)
      G = reverse!(G)
      # Update norms/scales from the updated Gram matrices:
      sL, aL, uL, wL = _scales_and_norms(G, p, uni)
      # Dualize again
      # JJ[steps[k][l]] *:= uniformizer^(-sL[k]);
      for k in 1:length(steps)
        for l in 1:length(steps[k])
          for u in 1:ncols(JJ)
            JJ[steps[k][l], u] = JJ[steps[k][l], u] * uni^(-sL[k])
          end
        end
      end

      # Update Gram matrices
      for k in 1:length(G)
        B = sub(JJ, steps[k][1]:steps[k][length(steps[k])], 1:ncols(JJ))
        G[k] = B * gram_matrix(ambient_space(L)) * B'
        @assert nrows(G[k]) in [1, 2]
      end
    end
  end

  @assert all(k -> nrows(G[k]) in [1,2], 1:length(G))
  return G, JJ
end

# Returns true iff BONG is a good BONG at p, as defined by C. Beli.
function isgood_bong(bong, p)
  v = all(valuation(bong[i], p) <= valuation(bong[i + 2], p) for i in 1:(length(bong)- 2))
  return v
end

function has_propertyA(L, p)
  @assert isdyadic(p)
  rL, sL, wL, aL = _genus_symbol_kirschmer(L, p).data
  nL = [valuation(aL[i], p) for i in 1:length(aL)]
  r = maximum(rL)
  if r > 2
    @vprint(:GenRep,1,"""Property A is violated over dyadic prime:
  There is a $(r)-dimensional Jordan component\n""")
  end

  # genus: rL, sL, wL, aL note that aL only contains valuations
  for i in 1:length(sL)
    for j in (i + 1):length(sL)
      if !((0 < nL[j] - nL[i]) && (nL[j] - nL[i] < 2*(sL[j] - sL[i])))
        @vprint(:GenRep, 1, """Property A is violated over dyadic prime:
  Violated at $(i) $(j) (norm/scale valuations do not fit)\n""")
        return false
      end
    end
  end
  return true
end

function G_function(a, V, g, p)
# Use LocalMultiplicativeGroupModSquares to calculate in F*/(O^*)^2
# (or F^*/(F^*)^2 -- the last component of the result is the p-valuation
# mod 2, and  F^*/(F^*)^2 =  F*/(O^*)^2 \times C_2.)
#
# Also we expect V, g = local_multiplicative_group_modulo_scares(p)
#
# cf. Beli 2003, Def. 4.
  O = order(p)
  K = nf(O)
  e = valuation(order(p)(2), p)
  R = valuation(a, p)
  d = relative_quadratic_defect(-a, p)

  if !is_in_A(a, p)
    @vprint :GenRep 2 "G_function case F\n"
    return N_function(-a, g, p)
  elseif g\(K(-1//4)) == g\(a)
    @vprint :GenRep 2 "G_function case G\n"
    return sub(V, gens(V)[1:ngens(V) - 1])
  elseif valuation(-4 * a, p) == 0 && islocal_square(-4 * a, p)
    @vprint :GenRep 2 "G_function case G\n"
    return sub(V, gens(V)[1:ngens(V) - 1])
  elseif R > 4 * e
    @vprint :GenRep 2 "G_function case H\n"
    return sub(V, [g\(a)])
  elseif 2*e < R && R <= 4 * e
    if d <= 2 * e - R//2
      @vprint :GenRep 2 "G_function case A\n"
      O = _one_plus_power_of_p(R + d - 2*e, V, g, p)
      return _intersect(N_function(-a, g, p), _sum(O, sub(V, [g\(a)])))
    else
      @vprint :GenRep 2 "G_function case B\n"
      @assert R % 2 == 0
      W, mW = _one_plus_power_of_p(div(R, 2), V, g, p)
      #@assert length(rels(W)) == 0
      return sub(V, append!([g\(a)], mW(w) for w in gens(W)))
    end
  elseif R <= 2 * e
    if d  <= e - R//2
      @vprint :GenRep 2 "G_function case C\n"
      return N_function(-a, g, p)
    elseif (e - R//2 < d) && (d <= 3 * e//2 - R//4)
      @assert R % 2 == 0
      @vprint :GenRep 2 "G_function case D\n"
      return _intersect(N_function(-a, g, p),
                        _one_plus_power_of_p(div(R, 2) + d - e, V, g, p))
    else
      @vprint :GenRep 2 "G_function case E\n"
      # Attention! Use the floor function wherever Beli writes stuff in square
      # brackets. This is due to his citing Hsia's papers, which have this
      # convention.
      return _one_plus_power_of_p(e - floor(Int, e//2 - R//4), V, g, p)
    end
  else
    throw(error("This should never happen"))
  end
end

# Returns the valuation of the relative quadratic defect of a at the prime
# ideal p, as defined by C. Beli
function relative_quadratic_defect(a, p)
  q = quadratic_defect(a, p)
  if q isa PosInf
    return inf
  else
    return q - valuation(a, p)
  end
end

# cf. Lemma 3.5 in Beli: A is the set of all a \in F^* for which either
# (-a\not\in (F^*)^2 and D(-a) \subseteq O_p) or (-a\in (F^*)^2 and ord(a)\geq
# -2e).
function is_in_A(a, p)
  e = valuation(order(p)(2), p)
  q = quadratic_defect(-a, p)
  if q isa Int
    if q >= 0
      return true
    end
  end

  if q isa PosInf
    if valuation(a, p) >= -2 * e
      return true
    end
  end

  return false
end

function N_function(a, g, p)
  #  g, ginv is the mapping obtained by local_multiplicative_group_modulo_squares(p).
  b = g\(a)

  V = parent(b)

  q = quadratic_defect(a, p)

  # project a into F^*/(O^*)^2:

  if q isa PosInf
    # treat the squares separately:
    return sub(V, gens(V))
  end

  B = gens(V)

  # cf. paragraph before 1.2 in Beli 2003:
  # N(a) := N(F(\sqrt(a))/F), i.e. the values of the norm mapping of the quadratic extension
  # and N(a) is the orthogonal complement of (<a>(F^*)^2).
  preimages = Vector{Int}(undef, length(B))
  for i in 1:length(B)
    b = B[i]
    preim = g(b)
    # convert the Hasse symbol (as an element of the multiplicative group C_2)
    # to an element of the additive group F_2:
    preimages[i] = hilbert_symbol(a, preim, p) == 1 ? 0 : 1
  end
  kernel_gen = elem_type(V)[ B[i] for i in 1:length(B) if preimages[i] == 0 ]
  # this fails if a is a square
  j = minimum(i for i in 1:length(B) if preimages[i] == 1)
  for i in 1:length(B)
    if preimages[i] == 1 && i != j
      push!(kernel_gen, B[i] + B[j])
    end
  end

  S, mS = sub(V, kernel_gen)
  #@assert length(rels(S)) == 0 # free
  return S, mS
end

# Map an idele into the ray class group associated to L. The parameter idele
# must be a list of tuples <p, a_p>, where p is a prime of base_ring(L), and a_p
# is an element of K^* (interpreted as an element of the completion K^*_p). The
# parameter atinfinity can be a list of tuples <v, +1 or -1>, where v is an
# element of real_places(nf(base_ring(L))). All places, finite or infinite, which
# are unspecified are interpreted as 1.}
function _map_idele_into_class_group(mRCG, idele, atinfinity::Vector{Tuple{InfPlc, Int}} = Tuple{InfPlc, Int}[])
  R = order(base_ring(codomain(mRCG)))
  F = nf(R)
  IP = defining_modulus(mRCG)[2]
  the_idele_inf = [1 for i in IP]
  if length(atinfinity) != 0
    for pl in atinfinity
      if pl[1] in IP
        the_idele_inf = pl[2]
      end
    end
  end
  M = defining_modulus(mRCG)[1]

  # The finite places we need to make congruent to 1 in order to be able to map into the class group:
  rayprimes = collect(keys(mRCG.fact_mod))
  exponents = Int[mRCG.fact_mod[p] for p in rayprimes]
  factors = ideal_type(R)[rayprimes[i]^exponents[i] for i in 1:length(rayprimes)]
  the_idele = [ one(F) for p in rayprimes ]
  for i in idele
    j = findfirst(isequal(i[1]), rayprimes)
    if j isa Int # found
      the_idele[j] = i[2]
    else
      #throw(error("Impossible?"))
      # ignore this
    end
  end
  # So, only the ideals D that are coprime to the ray lie in the codomain of mRCG and can be mapped
  #
  # The ideles we are considering here are considered to be representatives of
  # classes of ideles (modulo (F^*_S)*J^L, where F^*_S is the set of principal
  # ideles which are totally positive at the infinite places where the
  # envelopping space V of L is anisotropic, and where J^L is the set of ideles
  # which lie in the spinor norm of L at all finite places.  So, we may modify
  # the ideles in the following two ways without changing their class:
  # 1 multiply every component by an element of F^*
  # 2 multiply any individual component (at p) by an element of the spinor
  #   norms of L_p. In particular, any component is only relevant modulo
  #   squares.
  #
  # Thus we can achieve that the idele in question is congruent to 1 modulo the
  # ray M that defines the RCG
  #
  # The idele is then interpreted as the ideal \prod_p{p^Valuation(idele_p, p)}, which can be mapped into the class group.
  #
  # first, we need to be able to invert modulo the divisor of RCG
  #
  if any(valuation(the_idele[i], rayprimes[i]) != 0 for i in 1:length(the_idele))
    # We need to correct with an element of F^* which is positive at the relevant infinite places.
  # We will first correct only the finite places ...
    s = approximate([-valuation(the_idele[i], rayprimes[i]) for i in 1:length(the_idele)], rayprimes)
  else
    s = one(F)
  end

  # Now s * Idele is a unit at the RayPrimes
  # Now rescale once more with some t in F^* such that t*s*Idele = 1 mod MM and has the correct signs.

  sidele = Vector{elem_type(R)}(undef, length(rayprimes))
  for i in 1:length(the_idele)
    Q, mQ = quo(R, rayprimes[i]^exponents[i])
    sidele[i] = mQ\mQ(FacElem(s * the_idele[i]))
    # The quotient map can be applied on factored elements
  end

  if length(rayprimes) == 0
    x = one(R)
  else
    x = crt(sidele, factors)
    Q, mQ = quo(R, M)
    x = mQ\inv(mQ(x)) # x = invmod(x, M)
  end

  sgns = [ sign(s, IP[j]) * the_idele_inf[j] for j in 1:length(IP)]

  A, _exp, _log = infinite_primes_map(R, IP, M)
  t = x * (1 + _exp(A([ sgns[j] == sign(x, IP[j]) ? 0 : 1 for j in 1:length(IP)])))
  @assert x - t in M
  @assert all(sign(t, IP[j]) == sgns[j] for j in 1:length(IP))
  #t = crt(M, IP, x, sgns)

  s = s * t

  # Check if everything is ok.
  @hassert :GenRep 1 all(isone(quo(R, factors[k])[2](FacElem(s * the_idele[k]))) for k in 1:length(the_idele))
  @hassert :GenRep 1 all(sign(s * the_idele_inf[j], IP[j]) == 1 for j in 1:length(IP))

  # We first interpret it as the ideal which will actually have to be mapped:
  # i.e., we just collect the p-valuations at the noncritical places (p notin RayPrimes):

  _temp1 = [ idele[j][1] for j in 1:length(idele)]
  _temp2 = fmpz[ valuation(idele[j][2], idele[j][1]) for j in 1:length(idele)]
  ide = s * R
  _to_map = FacElem(Dict(numerator(ide) => fmpz(1), (denominator(ide) * R) => fmpz(-1)))
  if length(_temp1) != 0
    _to_map = _to_map * FacElem(_temp1, _temp2)
  end

  _to_map = FacElem(factor_coprime(_to_map))

  return g = mRCG\_to_map
end

function _compute_ray_class_group(L)
  R = base_ring(L)
  F = number_field(R)

  Gens = Tuple{ideal_type(R), elem_type(F)}[]

  rayprimes = ideal_type(R)[]
  M = 1 * R
  MM = ideal_type(R)[]
  Mfact = Dict{ideal_type(R), Int}()

  for p in bad_primes(L, even = true)
    spinors, V, g, exactlytheunits = spinor_norm(L, p)
    # we only need to carry around those finite places where the Spinor norm is
    # not exactly the units:
    if !exactlytheunits
      @vprint :GenRep 2 """Found a prime over $(minimum(p)) where the spinor
                           norm is not exactly the units of the order.
                             dim(spinors)=$(length(spinors)),
                             dim(LocalMultGrpModSq)=$(ngens(V))"""
      push!(rayprimes, p)
      # A basis of the spinor norm of L at p, when viewed (modulo squares) as an F_2-vector space
      b = spinors
      Gens = append!(Gens, [(p, g(a)) for a in b])
      @assert all(valuation(g(a), p) in [0, 1] for a in b)
      # We rely on local_multiplicative_group_modulo_squares (g) mapping to 0-
      # and 1-valued elements, and on the last base vector of V to represent a
      # uniformizer of p.
      I = p^(1 + 2 * valuation(R(2), p))
      M = M * I
      Mfact[p] = 1 + 2 * valuation(R(2), p)
      push!(MM, I)
    end
  end
  @assert issetequal(support(M), rayprimes)

  # Now M contains a ray M and MM is the support of this ray.
  # We now compute the indefinite real places of L
  inf_plc = [v for v in real_places(F) if !isisotropic(L, v)]
  # Now get the ray class group of M, inf_plc.
  return ray_class_group(M, inf_plc, lp = Mfact)..., Gens
end

function _get_critical_primes(L, mRCG, inf_plc, mQ, full = true)
  # Set full = true to calculate one prime for every element of
  # domain(mQ), instead of just one prime for every generator.
  # This is needed for indefinite lattices.

  R = base_ring(L)
  F = nf(R)
  critical_primes = ideal_type(R)[]
  Q = codomain(mQ)

  bad = prod(bad_primes(L, even = true))

  maxnorm = 50
  goodprimes = [ p for p in prime_ideals_up_to(R, maxnorm) if iscoprime(p, bad)]
  p_ind = 1

  grp_elements_to_primes = Dict{ideal_type(R), elem_type(Q)}()

  if full
    cur_list = elem_type(Q)[id(Q)]
    while length(cur_list) < length(Q)
      while p_ind > length(goodprimes)
        maxnorm = floor(Int, maxnorm * 1.2)
        goodprimes = ideal_type(R)[ p for p in prime_ideals_up_to(R, maxnorm) if iscoprime(p, bad)]
      end
      p = goodprimes[p_ind]
      h = mQ(_map_idele_into_class_group(mRCG, [(p, F(uniformizer(p)))]))
      if !haskey(grp_elements_to_primes, h)
        grp_elements_to_primes[p]= h
      end
      oldsize = length(cur_list)
      if !(h in cur_list)
        push!(cur_list, h)
      end
      push!(critical_primes, p)
      p_ind = p_ind + 1
    end
  else
    throw(NotImplemented())
  end

  return critical_primes
end

function _spinor_generators(L, C, mod_out = elem_type(codomain(C.mQ))[])
  R = base_ring(L)
  bad = ideal_type(R)[ p for p in bad_primes(L) if !ismodular(L, p)[1] ]
  S, mS = sub(codomain(C.mQ), mod_out)
  n = order(codomain(C.mQ))
  gens = ideal_type(R)[]
  tmp_gens = copy(mod_out)
  R = base_ring(L)
  p = 2
  while order(S) != n
    p = next_prime(p)
    lp = ideal_type(R)[ p[1] for p in prime_decomposition(R, p) if !(p[1] in bad) ]
    for P in lp
      g = _map_idele_into_class_group(C.mR, [(P, uniformizer(P))])
      h = C.mQ(g)
      if !haspreimage(mS, h)[1]
        push!(gens, P)
        push!(tmp_gens, h)
        S, mS = sub(codomain(C.mQ), tmp_gens)
      end
      if order(S) == n
        break
      end
    end
  end
  return gens
end

# TODO: Enable use_auto
function neighbours(L::QuadLat, p; call = stdcallback, use_auto = true, max = inf)
  R = base_ring(L)
  F = nf(R)
  @req R == order(p) "Incompatible arguments"
  @req isprime(p) "Ideal must be prime"
  ok, rescale = ismodular(L, p)
  @req ok "The lattice must be locally modular"
  @req isisotropic(L, p) "The lattice must be locally isotropic"
  e = valuation(R(2), p)
  @req e == 0 || valuation(norm(L), p) >= e "The lattice must be even"
  B = local_basis_matrix(L, p, type = :submodule)
  n = nrows(B)
  if F isa AnticNumberField
    @assert nbits(minimum(p)) < 60
    k, h = ResidueFieldSmall(R, p)
  else
    k, h = ResidueField(R, p)
  end
  hext = extend(h, F)
  pi = uniformizer(p)
  piinv = anti_uniformizer(p)
  form = gram_matrix(ambient_space(L), B)
  if rescale != 0
    form = piinv^rescale * form
  end
  pform = map_entries(hext, form)

  if use_auto
    G = automorphism_group_generators(L)
    @hassert :GenRep 1 all(g -> g * gram_matrix(ambient_space(L)) * transpose(g) == gram_matrix(ambient_space(L)), G)
    Binv = inv(B)
    adjust_gens = eltype(G)[B * g * Binv for g in G]
    @hassert :GenRep 1 all(g -> g * form * transpose(g) == form, adjust_gens)
    adjust_gens_mod_p = dense_matrix_type(k)[map_entries(hext, g) for g in adjust_gens]
    adjust_gens_mod_p = dense_matrix_type(k)[x for x in adjust_gens_mod_p if !isdiagonal(x)]
    @hassert :GenRep 1 all(g -> g * pform * transpose(g) == pform, adjust_gens_mod_p)
    q = order(k)
    if length(adjust_gens_mod_p) > 0
      _LO = line_orbits(adjust_gens_mod_p)
      LO = Vector{eltype(k)}[x[1] for x in _LO]
      @vprint :GenRep 2 "Checking $(length(LO)) representatives (instead of $(div(order(k)^n - 1, order(k) - 1)))\n"
    else
      @vprint :GenRep 2 "Enumerating lines over $k of length $n\n"
      LO = enumerate_lines(k, n)
    end
  else
    @vprint :GenRep 2 "Enumerating lines over $k of length $n\n"
    LO = enumerate_lines(k, n)
  end

  result = typeof(L)[]

  pMmat = _module_scale_ideal(p, pseudo_matrix(L))

  # TODO: This is too slow
  _dotk(u, v) = (matrix(k, 1, n, u) * pform * matrix(k, n, 1, v))[1, 1]

  _dotF(u, v) = (matrix(F, 1, n, u) * form * matrix(F, n, 1, v))[1, 1]

  keep = true
  cont = true
  found = false

  for w in LO
    dotww = _dotk(w, w)
    if dotww != 0
      continue
    end
    if iszero(matrix(k, 1, n, w) * pform)
      continue
    end

    x = [ hext\e for e in w]

    nrm = _dotF(x, x)
    val = valuation(nrm, p)
    @assert val > 0
    if val <= e
      continue # can only happen if p is even
    elseif val == e + 1
      # make val > e + 1
      bas = [zero(k) for i in 1:n]
      r = 0
      for i in 1:n
        bas[i] = one(k)
        if _dotk(w, bas) != 0
          r = i
          break
        end
        bas[i] = zero(k)
      end
      @assert r != 0
      u = [zero(F) for i in 1:n]
      u[r] = one(F)
      a = (h\(hext((nrm//(2 * pi * _dotF(x, u))))))
      @. x = x - a * pi * u
      @assert valuation(_dotF(x, x), p) >= e + 2
    end
    found = true

    # normalize x

    kk = 0
    for i in 1:n
      if hext(x[i]) != 0
        kk = i
        break
      end
    end
    @assert kk != 0
    x = Ref(h\(inv(hext(x[kk])))) .* x
    if e != 0
      x = Ref(1 + h\(hext((x[kk] - 1)//pi)) * pi) .* x
      @assert valuation(x[kk] - 1, p) >= 2
    end

    xF = matrix(F, 1, n, x) * form
    mm = 0
    for r in 1:n
      if r != kk && !iszero(xF[1, r]) && valuation(xF[1, r], p) == 0
        mm = r
        break
      end
    end
    @assert mm != 0

    _U = [zero(F) for i in 1:n]
    _U[mm] = one(F)

    _tt = [ Ref(piinv) .* x, Ref(0 * pi) .* _U]
    for i in 1:n
      if i != kk && i != mm
        _t = [zero(F) for i in 1:n]
        _t[i] = one(F)
        push!(_tt, _t .- (Ref(h\hext(xF[1, i]//xF[mm])) .* _U))
      end
    end
    VV = matrix(F, n, n, reduce(vcat, _tt))
    V = VV * B
    LL = lattice(ambient_space(L), _sum_modules(pMmat, pseudo_matrix(V)))

    @hassert :GenRep 1 islocally_isometric(LL, L, p)

    if !(call isa Bool)
      keep, cont = call(result, LL)
    end
    if keep
      push!(result, LL)
    end

    if !cont || length(result) >= max
      break
    end
  end
  if !found
    @warn "L_p/pL_p has no suitable isotropic vector!"
  end
  return result
end

function iterated_neighbours(L::QuadLat, p; use_auto = true, max = inf, mass = -one(FlintQQ))
  @req isdefinite(L) "Lattice must be definite"
  result = typeof(L)[ L ]
  i = 1

  use_mass = mass > 0

  local found::fmpq

  if mass >= 0
    found = 1//automorphism_group_order(L)
  end

  while (i <= length(result)) && (length(result) < max) && (!use_mass || found < mass)
    # keep if not isometric, continue until the whole graph has been exhausted.
    callback = function(res, M)
      keep = all(LL -> !isisometric(LL, M)[1], vcat(res, result))
      return keep, true;
    end
    N = neighbours(result[i], p, call = callback, use_auto = use_auto, max = max - length(result))
    if use_mass && !isempty(N)
      found = found + sum(fmpq[1//automorphism_group_order(LL) for LL in N])
      perc = Printf.@sprintf("%2.1f", Float64(found//mass) * 100)
      @vprint :GenRep 1 "#Lattices: $(length(result)), Target mass: $mass. Found so far: $found ($perc%)\n"
    end
    append!(result, N)
    i = i + 1
  end
  return result
end

function _scales_and_norms(G, p, uni)
  R = order(p)
  K = nf(R)
  e = ramification_index(p)
  aL = elem_type(K)[]
  uL = Int[]
  wL = Int[]
  sL = Int[ minimum(Union{Int, PosInf}[iszero(x) ? inf : valuation(x, p) for x in g]) for g in G ]
  for i in 1:length(G)
    # similar, but not equal to, the code for obtaining a genus symbol
    # (for the genus symbol, we consider a scaling L^(s(L_i)) in O'Meara's notation)

    GG = G[i]
    D = diagonal(GG)
    if e + sL[i] <= minimum(Union{PosInf, Int}[iszero(d) ? inf : valuation(d, p) for d in D])
      push!(aL, elem_in_nf(uni^(e + sL[i])))
    else
      _, b = findmin([valuation(x, p) for x in D])
      push!(aL, D[b])
    end
    push!(uL, valuation(aL[i], p))
    @assert uL[i] == valuation(norm(quadratic_lattice(nf(R), identity_matrix(K, nrows(G[i])), gram_ambient_space = GG)), p)
    push!(wL, min(e + sL[i], minimum(uL[i] + quadratic_defect(d//aL[i], p) for d in D)))
  end
  return sL, aL, uL, wL
  # scales, norm generators, norm valuations, weight valuations of a (Jordan) decomposition of L
end

function __ismaximal_norm_splitting(gram_matrices, scales, norms, p)
  #  Scales: list of valuations of scales
  #  Norms: list of generators of norms
  #  occurring in a Genus symbol of L at p, as calculated by GenusSymbol(L, p).
  normsval = Int[valuation(norms[i], p) for i in 1:length(norms) ]
  # test if each component is either unary or binary:
  k = findfirst(i -> !(ncols(gram_matrices[i]) in [1, 2]), 1:length(gram_matrices))
  if k !== nothing
    @vprint :GenRep 2 "not maximal norm splitting: components are not all unary or binary\n";
    return false, -k, []
  end
  # test if binary components are modular:
  for i in 1:length(gram_matrices)
    if ncols(gram_matrices[i]) != 1 && valuation(det(gram_matrices[i]), p) != 2*scales[i]
      @vprint :GenRep 2 "not maximal norm splitting: at least one of the binary components is not modular\n"
      return false, -i, []
    end
  end
  # test if sL[1] \supseteq sL[2] \supseteq ... \supseteq sL[#sL]:
  for i in 1:length(scales) - 1
    if scales[i] > scales[i + 1]
      throw(Error("Your lattice is weird"))
      return fail, 0, []
    end
  end

  NU, _ = _norm_upscaled(gram_matrices, p)
  # test if nL[i] = n(L^{sL[i]}):
  fail = []
  for i in 1:length(scales)
    @assert NU[i] <= normsval[i]
    if NU[i] < normsval[i]
      @vprint :GenRep 2 "not maximal norm splitting: norm condition at $i\n"
      push!(fail, i)
    end
  end

  if length(fail) > 0
    return false, fail[1], fail
  end
  return true, 0, []
end

function _ismaximal_norm_splitting(G, p)
  sL, aL, _, _ = scales_and_norms(G, p, uniformizer(p))
  return __is_maximal_norm_splitting(G, sL, aL, p)
end

function _norm_upscaled(G, p)
  # calculate Norm(L^{scale(L_i)}), for all i.
  # for the definition of L^{a}, cf. § 82 I, p.237, in O'Meara
  #
  # G is a splitting of a LatticeModule L into Gram matrices,
  # which can e.g. be calculated by JordanDecomposition(), but need not
  # be a Jordan decomposition.
  t = length(G)
  sL = Int[ minimum(Union{Int, PosInf}[iszero(g[i, j]) ? inf : valuation(g[i, j], p) for j in 1:ncols(g) for i in 1:j]) for g in G]
  e = ramification_index(p)
  uni = elem_in_nf(uniformizer(p))
  aL = []
  uL = []
  wL = []
  for i in 1:t
    GG = diagonal_matrix([ j < i ? uni^(2*(sL[i] - sL[j])) * G[j] : G[j] for j in 1:t])
    # the norm is 2*Scale + <ideals generated by the diagonals>, cf. § 94 O'Meara.
    D = diagonal(GG)
    # Is 2*s(L_i) eq n_i? (We always have scale \supseteq norm \supseteq 2*scale)
    # To find norm generator, pick a generator of 2*Scale or <diagonal elements>
    if e + sL[i] <= minimum(valuation(d, p) for d in D)
      # 2*scale is the bigger ideal
      push!(aL, uni^(e + sL[i]))
    else
      # diagonal elements for the bigger ideal
      _, b = findmin([valuation(x, p) for x in D])
      push!(aL, D[b])
    end
    push!(uL, valuation(aL[i], p))
  end
 # uL: the valuations of the upscaled norms
 # aL: the generators of the upscaled norms
 return uL, aL
end

function _make_bong_dim_2(L, p)
  # cf. Beli, Lemma 3.3
  # return a BONG of a 2-dimensional lattice.

  ismod, r = ismodular(L, p)
  @assert ismod && rank(L) == 2
  pi = uniformizer(p)
  R = order(p)
  K = nf(R)
  A = gram_matrix(ambient_space(L))
  # is the (1,1) entry a norm generator?
  if valuation(A[1, 1], p) != valuation(norm(L), p)
    # swap stuff around so that the (1,1) entry will be a norm generator
    b = 0
    for i in 1:rank(A)
      if valuation(A[i, i], p) == valuation(norm(L), p)
        b = i
        break
      end
    end

    if b != 0
      S = local_basis_matrix(L, p)
      swap_rows!(S, 1, 2)
    else
      S = local_basis_matrix(L, p)
      for j in 1:ncols(S)
        S[1, j] = S[1, j] + S[2, j]
      end
    end
    A = S * gram_matrix(ambient_space(L)) * transpose(S)
    Lold = L
    L = quadratic_lattice(K, identity_matrix(K, nrows(A)), gram_ambient_space = A)
    @assert valuation(A[1, 1], p) == valuation(norm(L), p)
  end
  n = A[1, 1]
  # Some final checks:
  disc = A[1, 1] * A[2, 2] - A[2, 1] * A[1, 2]
  val = valuation(disc, p)
  @assert iseven(val)
  @assert valuation(disc, p) == 2*r
  bong = [n, disc//n]
end

function _one_plus_power_of_p(k, V, g, p)
  # See Beli 2003, Def. 1.
  # We expect V, g = local_multiplicative_group_modulo_squares(p)
  r = ngens(V)
  it = Iterators.product([collect(0:1) for i in 1:r]...)
  S = [ g(V(collect(v))) for v in it ]
  SS = [ s for s in S if relative_quadratic_defect(s, p) >= k ]
  return sub(V, [g\(s) for s in SS])
end

function _intersect(_V, _W)
  V, mV = _V
  W, mW = _W
  return intersect(mV, mW)
end

function _sum(_V, _W)
  V, mV = _V
  W, mW = _W
  bigspace = codomain(mV)
  return sub(bigspace, append!([mV(g) for g in gens(V)], (mW(g) for g in gens(W))))
end

function beli_correction!(L, G, JJ, steps, i, j, p)
  # this is helper procedure for GoodBONG().
  # Correct the Jordan components with number i and j of L from the Jordan
  # decomposition given in G and JJ, in such a way that the new component with
  # number i has maximal norm.
  @assert length(steps[i]) == 2
  NU = _norm_upscaled(G, p)
  if nrows(G[j]) == 2
    # assert orthogonality of the vectors in JJ[Steps[i]] and those in
    # JJ[Steps[j]], i.e. those that make up G[i] and G[j]:
    #B = sub(JJ, steps[i][1]:(steps[i][1] + steps[i][2] - 1), 1:ncols(JJ))
    B = zero_matrix(base_ring(JJ), 2, ncols(JJ))
    for k in 1:ncols(JJ)
      B[1, k] = JJ[steps[i][1], k]
      B[2, k] = JJ[steps[i][2], k]
    end
    C = sub(JJ, steps[j][2]:steps[j][2], 1:ncols(JJ))
    @assert iszero(C * gram_matrix(ambient_space(L)) * B')
  end

  K = nf(base_ring(L))

  # assert NU[i] lt Valuation(Norm(LatticeModule(G[i])), p); // otherwise, the norm condition is violated
  # assert NU[j] le Valuation(Norm(LatticeModule(G[j])), p); // "<=" must always hold, otherwise something about the lattice is broken
#  */
#
  #  We will need the original Gram matrix for re-orthogonalization:
  GI = inv(G[i])
  @assert length(steps[j]) in [1, 2]
  # if G[j] is two-dimensional, make sure that a norm generator is on the [1,1]
  # position:
  if ncols(G[j]) == 2 && valuation(G[j][1, 1], p) > valuation(G[j][2, 2], p)
    temp = JJ[steps[j][1]]
    swap_rows!(JJ, steps[j][1], steps[j][2])
    B = zero_matrix(base_ring(JJ), 2, ncols(JJ))
    for u in 1:ncols(JJ)
      B[1, u] = JJ[steps[j][1], u]
      B[2, u] = JJ[steps[j][2], u]
    end
    G[j] = B * gram_matrix(ambient_space(L)) * B'
  end

  for u in 1:ncols(JJ)
    JJ[steps[i][1], u] = JJ[steps[i][1], u] + JJ[steps[j][1], u]
  end
  # update Gram matrix for component #i:
  B = zero_matrix(base_ring(JJ), 2, ncols(JJ))
  for u in 1:ncols(JJ)
    B[1, u] = JJ[steps[i][1], u]
    B[2, u] = JJ[steps[i][2], u]
  end

  G[i] = B * gram_matrix(ambient_space(L)) * B'
  x = sub(JJ, steps[i][1]:steps[i][1], 1:ncols(JJ))
  y = sub(JJ, steps[i][2]:steps[i][2], 1:ncols(JJ))

  # Ansatz: v = JJ[Steps[j][1]] + lambda*x + mu*y

  # re-orthogonalize the first basis vector of G[j]:
  v = matrix(K, 1, 2, [-G[j][1, 1], 0]) * GI
  # G[j][1,1] is always non-zero (the lattice is definite)
  # assert integrality at p:
  @assert all(k -> valuation(v[k], p) >= 0, 1:ncols(v))

  # JJ[steps[j][1]] +:= v[1]*JJ[steps[i][1]] + v[2]*JJ[steps[i][2]];
  for u in 1:ncols(JJ)
    JJ[steps[j][1], u] = JJ[steps[j][1], u] + v[1] * JJ[steps[i][1],u] + v[2] * JJ[steps[i][2], u]
  end

  # if applicable, re-orthogonalize the second basis vector of G[j]:
  if ncols(G[j]) == 2
    w = matrix(K, 1, 2, [-G[j][1, 2], 0]) * GI
    # G[j][1,2] is always non-zero (or else the lattice would have been further
    # decomposed into two 1x1-lattices here)

    @assert all(k -> valuation(w[k], p) >= 0, 1:ncols(w))
    for u in 1:ncols(JJ)
      JJ[steps[j][1], u] = JJ[steps[j][1], u] + w[1] * JJ[steps[i][1],u] + w[2] * JJ[steps[i][2], u]
    end
  end

  B = sub(JJ, steps[j][1]:steps[j][length(steps[j])], 1:ncols(JJ))
  G[j] = B * gram_matrix(ambient_space(L)) * B'
end

################################################################################
#
#  Local multiplicative modulo squares
#
################################################################################

# Helper function to find y such that xy^2 is integral and Valuation(xy^2, p) = 0
function _square_rep_nice(x, p, piinv)
  x = x * denominator(x, order(p))^2
  v = valuation(x, p)
  @assert v >= 0 && iseven(v)
  if !iszero(v)
    x = x * piinv^v
  end
  return order(p)(x)
end

# Move this to a proper place
#
# TODO: Cache this in the dyadic case (on the lattice or the field)
mutable struct LocMultGrpModSquMap <: Map{GrpAbFinGen, GrpAbFinGen, HeckeMap, LocMultGrpModSquMap}
  domain::GrpAbFinGen
  codomain::AnticNumberField
  isdyadic::Bool
  p::NfAbsOrdIdl{AnticNumberField, nf_elem}
  e::nf_elem
  pi::nf_elem
  piinv::nf_elem
  hext::NfToFinFldMor{FqFiniteField}
  h::AbsOrdQuoMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}
  g::GrpAbFinGenToAbsOrdQuoRingMultMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}
  i::GrpAbFinGenMap
  mS::GrpAbFinGenMap

  function LocMultGrpModSquMap(K::AnticNumberField, p::NfAbsOrdIdl{AnticNumberField, nf_elem})
    R = order(p)
    @assert nf(R) === K
    @assert isabsolute(K)
    z = new()
    z.codomain = K
    z.p = p
    z.isdyadic = isdyadic(p)

    if !isdyadic(p)
      pi = elem_in_nf(uniformizer(p))
      k, h = ResidueField(R, p)
      hext = extend(h, K)
      e = elem_in_nf(h\non_square(k))
      G = abelian_group([2, 2])

      z.domain = G
      z.e = e
      z.pi = pi
      z.hext = hext
      return z
    else
      pi = elem_in_nf(uniformizer(p))
      e = ramification_index(p)
      dim = valuation(norm(p), 2) * e + 2
      #V = vector_space(F, dim)
      I = p^(2*e + 1)
      Q, h = quo(R, I)
      U, g = unit_group(Q)
      M, i = quo(U, 2, false)
      SS, mSS = snf(M)
      @assert SS.snf == fmpz[2 for i in 1:(dim - 1)]
      #@assert ngens(S) == dim - 1
      piinv = anti_uniformizer(p)
      G = abelian_group([2 for i in 1:dim])
      z.domain = G
      z.pi = pi
      z.piinv = piinv
      z.h = h
      z.g = g
      z.i = i
      z.mS = mSS
      return z
    end
  end
end

domain(f::LocMultGrpModSquMap) = f.domain

codomain(f::LocMultGrpModSquMap) = f.codomain

function show(io::IO, f::LocMultGrpModSquMap)
  print(io, "Map for local unit group modulo squares at\n")
  print(IOContext(io, :compact => true), f.p, "\n")
  print(io, "from\n")
  print(IOContext(io, :compact => true), "  ", f.domain, "\n")
  print(io, "to\n")
  print(IOContext(io, :compact => true), "  ", f.codomain)
end

function image(f::LocMultGrpModSquMap, x::GrpAbFinGenElem)
  @assert parent(x) == f.domain
  K = f.codomain
  if !f.isdyadic
    if iszero(x.coeff[1])
      y = one(K)
    else
      y = f.e
    end
    if isone(x.coeff[2])
      y = y * f.pi
    end
    return y
  else
    S = codomain(f.mS)
    G = domain(f)
    y = elem_in_nf(f.h\(f.g(f.i\(f.mS(S(fmpz[x.coeff[i] for i in 1:(ngens(G) - 1)]))))))
    if x.coeff[ngens(G)] == 1
      y = y * f.pi
    end
    return y
  end
end

function preimage(f::LocMultGrpModSquMap, y::nf_elem)
  @assert parent(y) == f.codomain
  if !f.isdyadic
    v = valuation(y, f.p)
    if issquare(f.hext(y//f.pi^v))[1]
      fir = 0
    else
      fir = 1
    end
    return f.domain([fir, v])
  else
    v = valuation(y, f.p)
    w = f.mS\(f.i(f.g\(f.h(_square_rep_nice(y * f.pi^v, f.p, f.piinv)))))
    G = domain(f)
    return G(push!(fmpz[w.coeff[i] for i in 1:(ngens(domain(f)) - 1)], v))
  end
end

## Given a prime ideal over some number field K, this returns a vector space
## over GF(2) isomorphic to K_p^/(K_p^*)^2 and a map representing the isomorphism
##
function local_multiplicative_group_modulo_squares(p)
  K = nf(order(p))
  f = LocMultGrpModSquMap(K, p)
  return domain(f), f
end

function non_square(F::Union{GaloisField, FqFiniteField})
  r = rand(F)
  while iszero(r) || issquare(r)[1]
    r = rand(F)
  end
  return r
end

################################################################################
#
#  Binary quadratic forms
#
################################################################################

function __colon_raw(K, a, b)
  d = degree(K)
  bb = b
  B = inv(basis_matrix(a)) #fmpq_mat(basis_mat_inv(a, copy = false))
  M = zero_matrix(FlintQQ, d^2, d)
  for i = 1:d
    N = representation_matrix(bb[i])*B
    for s = 1:d
      for t = 1:d
        M[t + (i - 1)*d, s] = N[s, t]
      end
    end
  end
  M = sub(hnf(FakeFmpqMat(M), :upperright), 1:d, 1:d)
  N = inv(transpose(M))
  return N
end

function _genus_representatives_binary_quadratic_definite(L::QuadLat; max = inf, use_auto = true, use_mass = true)
  # The internal functions wants 1 in Q(K \otimes L)
  # So we rescale a bit

  V = rational_span(L)
  D = diagonal(V)
  _, i = findmin(abs.(norm.(D)))
  d = D[i]
  # Do G -> d * G
  _L = rescale(L, d)
  lat = _genus_representatives_binary_quadratic_definite_helper(_L; max = max, use_auto = use_auto, use_mass = use_mass)
  G = genus(L)
  res = typeof(L)[]
  for M in lat
    Mre = rescale(M, inv(d))
    @test genus(Mre) == G
    push!(res, Mre)
  end
  return res
end

function _genus_representatives_binary_quadratic_definite_helper(L::QuadLat; max = inf, use_auto = true, use_mass = true)
  # The strategy is to pass to the discriminant field F (which is a CM field)
  # and use that KL \cong (F, Tr/2)
  # Then in (F, Tr/2) we use Kirschmer, Pfeuffer and Körrner.

  K = base_field(L)
  L = lattice(quadratic_space(base_field(L), gram_matrix_of_rational_span(L)),
                              pseudo_matrix(identity_matrix(K, 2), coefficient_ideals(L)))

  @assert isdefinite(L)
  @assert rank(L) == 2

  #V = rational_span(_L)
  ## 0. Scale to have 1 in Q(V)
  #@show V

  #D = diagonal(V)
  #_, i = findmin(abs.(norm.(D)))
  #d = D[i]

  ## so G -> G/d

  #L = lattice(quadratic_space(base_ring(V), 1//d * gram_matrix(ambient_space(_L))), pseudo_matrix(_L))

  V = rational_span(L)

  # 1. Find the isometry with (F, Tr/2)
  @vprint :GenRep 1 "Determining isometry with CM field ... \n"
  K = base_ring(V)
  d = discriminant(V)
  de = denominator(d)
  @assert !issquare(de * d)[1]
  Kt, t = PolynomialRing(K, "t", cached = false)
  F, z = number_field(t^2 - de^2 * d, "z", cached = false)
  # TODO: Use automorphism_group (once implemented for relative extensions)
  a1, a2 = automorphisms(F)
  if a1(z) == z
    sigma = a2
  else
    sigma = a1
  end

  B = basis(F)
  phi(x, y) = K((x * sigma(y) + y * sigma(x))//2)
  G = matrix(K, 2, 2, [phi(B[1], B[1]), phi(B[1], B[2]), phi(B[2], B[1]), phi(B[2], B[2])])
  W = quadratic_space(K, G)
  fl, T = isequivalent_with_isometry(V, W)
  # Note that this is an isometry of KL with W
  @assert fl
  Tinv = inv(T)

  # Compute the absolute field, since this is where we do most of the work.

  Fabs, FabstoF = absolute_simple_field(F)

  KtoFabs = restrict(inv(FabstoF), base_field(F))

  sigmaabs = hom(Fabs, Fabs, FabstoF\(sigma(FabstoF(gen(Fabs)))))

  # 2. Transpose L to lattice of F.

  pb = pseudo_basis(L)
  # KL has basis a[1] for a in pb
  # So we only need to map the Z-basis of the coefficient ideals
  image_of_first = FabstoF\(T[1, 1] * B[1] + T[1, 2] * B[2])
  image_of_second = FabstoF\(T[2, 1] * B[1] + T[2, 2] * B[2])

  #@show [phi(FabstoF(a), FabstoF(b)) for a in [image_of_first, image_of_second] for b in [image_of_first, image_of_second]]

  basisofLinFabs = elem_type(Fabs)[]

  B = absolute_basis(pb[1][2])
  for b in B
    push!(basisofLinFabs, FabstoF\(F(b)) * image_of_first)
  end

  B = absolute_basis(pb[2][2])
  for b in B
    push!(basisofLinFabs, FabstoF\(F(b)) * image_of_second)
  end

  # Let M = <basisofLinFabs>_Z, this is a lattice in Fabs and we need to compute
  # its (left/right) order. We use __colon_raw on the basis matrix for this.

  EFabs = equation_order(Fabs)

  scaled_basisofLinFabs = elem_type(EFabs)[]

  O = Order(Fabs, __colon_raw(Fabs, basisofLinFabs, basisofLinFabs))

  z = zero_matrix(K, degree(O), degree(F))
  for i in 1:degree(O)
    elem_to_mat_row!(z, i, FabstoF(elem_in_nf(basis(O)[i])))
  end
  pm = pseudo_hnf(pseudo_matrix(z), :lowerleft)
  i = 1
  while iszero_row(pm.matrix, i)
    i += 1
  end
  pm = sub(pm, i:nrows(pm), 1:ncols(pm))

  OinF = Order(F, pm)

  LinFabs = fractional_ideal(O, basisofLinFabs)

  # O is the ring of multipliers/right oder of M
  # Now we compute CL_F/{ambiguous ideals of O_F with respect to F/K}

  @vprint :GenRep 1 "Compute representatives modulo ambiguous ideal classes ... \n"
  OFabs = maximal_order(O)
  amb_ideals = _ambigous_ideals_quotient_basis(OFabs, F, KtoFabs, FabstoF)
  @assert all(sigmaabs(I) == I for I in amb_ideals)
  OK = maximal_order(K)
  CK, mCK = class_group(OK)
  ideals_from_K = ideal_type(OFabs)[ KtoFabs(mCK(c)) for c in gens(CK)]
  CF, mCF = class_group(OFabs)
  Q, mQ = quo(CF, append!(elem_type(CF)[mCF\I for I in amb_ideals], elem_type(CF)[mCF\I for I in ideals_from_K]))
  repr = ideal_type(OFabs)[mCF(mQ\q) for q in Q]

  # There is a formula for the number of ambiguous ideal classes, this checks this:
  # UK, mUK = unit_group(OK)
  # UFabs, mUFabs = unit_group(OFabs)
  #@show order(quo(UFabs, append!([ mUFabs\OFabs(KtoFabs(elem_in_nf(mUK(u)))) for u in gens(UK) ], [mUFabs\OFabs(torsion_units_generator(Fabs))]))[1])
  #@show order(CF)
  #@show order(CK)
  #@show length(support(discriminant(maximal_order(F))))

  # repr are representatives for CL_F/{ambiguous ideals of O_F with respect to F/K}
  # This is isomorphic to gen(OFabs)/cls+(OFabs) via a -> b = a/\sigma(a).
  # But I also want to compute x = (x_p)_p with b_p = x_p O_p and
  # x_p \in F_p^1
  #
  # This is done in _translate_ideal
  #
  # _intersect_lattice_down computes a preimage under gen(O) -> gen(CL_F)
  #
  # The map gen(O) -> gen(CL_F) is not injective, but we can compute a set containing
  # representatives of the kernel elements.

  OF = maximal_order(F)
  repr_in_OF = ideal_type(OF)[ reduce(+, (OF(FabstoF(elem_in_nf(b))) * OF for b in basis(repr[i])), init = OF(0) * OF) for i in 1:length(repr)]

  c = conductor(O, OFabs) * OFabs
  _xps = Vector{elem_type(F)}[]
  _ps = ideal_type(OK)[]
  @vprint :GenRep 1 "Computing local units ...\n"
  for (p, _) in factor(minimum(image(FabstoF, c)))
    lQ = prime_decomposition(OF, p)
    lQ = [ (preimage(FabstoF, Q), e) for (Q, e) in lQ ]
    cppart = reduce(*, (Q^valuation(c, Q) for (Q, _) in lQ), init = 1 * OFabs)
    Rq, mq = quo(OFabs, cppart)
    Mu, f = multiplicative_group(Rq)
    elts = elem_type(F)[ FabstoF(elem_in_nf(mq\(f(m)))) for m in Mu]
    if length(elts) > 1
      push!(_ps, p)
      push!(_xps, elts)
    end
  end

  kernel_rep = fractional_ideal_type(O)[]

  if length(_xps) > 0
    IT = Iterators.product(_xps...)
    @vprint :GenRep 1 "There are $(length(IT)) potential kernel elements.\n"

    for it in IT
      LL = _intersect_lattice_down(collect(it)::Vector{elem_type(F)}, _ps, OinF)
      LLinabs = fractional_ideal(O, [FabstoF\(b) for b in absolute_basis(LL)])
      push!(kernel_rep, LLinabs)
    end
  end

  almost_res = fractional_ideal_type(O)[]

  index_divs = support(fmpq(index(O, OFabs)))

  for i in 1:length(repr_in_OF)
    I = repr_in_OF[i]
    Iabs = repr[i]
    N, xps, ps = _translate_ideal(I, Iabs, FabstoF, sigma, sigmaabs)
    @assert norm(N) == 1 * OK
    M = _intersect_lattice_down(xps, ps, OinF)

    Minabs = fractional_ideal(O, [FabstoF\(b) for b in absolute_basis(M)])
    if length(kernel_rep) > 0
      for K in kernel_rep
        push!(almost_res, K * Minabs * LinFabs)
      end
    else
      push!(almost_res, Minabs * LinFabs)
    end
  end

  @vprint :GenRep 1 "Now sieving the $(length(almost_res)) many candidates ...\n"
  # At the end we have to translate this relative ideals again and then back to
  # V. The things we want have the underscore _, but for assertion purposes we
  # do both.

  gram_ambient_space = gram_matrix(V)

  V = ambient_space(L)

  G = genus(L)

  res = typeof(L)[]

  cur_mass = zero(fmpq)

  if use_mass
    _mass = mass(L)
    @vprint :GenRep 1 "Using mass, which is $(_mass)\n"
  else
    _mass = one(fmpq)
    @vprint :GenRep 1 "Not using mass\n"
  end

  for N in almost_res
    generators = []
    for b in basis(N)
      binF = FabstoF(b)
      push!(generators, [Tinv[1, 1] * coeff(binF, 0) + Tinv[2, 1] * coeff(binF, 1),
                         Tinv[1, 2] * coeff(binF, 0) + Tinv[2, 2] * coeff(binF, 1)])
    end

    z = zero_matrix(K, length(generators), dim(V))
    for i in 1:length(generators)
      for j in 1:ncols(gram_ambient_space)
        z[i, j] = generators[i][j]
      end
    end
    pm = pseudo_hnf(pseudo_matrix(z), :lowerleft)
    i = 1
    while iszero_row(pm.matrix, i)
      i += 1
    end
    pm = sub(pm, i:nrows(pm), 1:ncols(pm))

    _pm = pseudo_matrix(matrix(pm), coefficient_ideals(pm))

    _new_cand = lattice(V, pm)

    if genus(_new_cand) != G
      continue
    end

    if any(T -> isisometric(T, _new_cand)[1], res)
      continue
    else
      push!(res, _new_cand)
      if length(res) >= max
        break
      end
      if use_mass
        cur_mass += fmpq(1, automorphism_group_order(_new_cand))
        if cur_mass > _mass
          error("Something very wrong in genus_representatives (mass mismatch")
        end
        if cur_mass == _mass
          break
        end
      end
    end
  end

  if max < inf && use_mass
    @assert _mass == cur_mass
  end

  return res
end

function _translate_ideal(I, Iabs, FtoFabs, sigma, sigmaabs)
  O = order(I)
  crit_primes = support(minimum(I))
  F = nf(O)
  xps = elem_type(F)[]
  sigmaI = sum(ideal(order(I), order(I)(FtoFabs(elem_in_nf(x)))) for x in basis(sigmaabs(Iabs)))
  M = I * inv(sigmaI)
  for p in crit_primes
    lQ = prime_decomposition(O, p)
    lQ_ = ideal_type(O)[ Q for (Q, _) in lQ]
    pvals = Int[]
    for Q in lQ_
      push!(pvals, valuation(I, Q))
    end

    zp = approximate(pvals, lQ_)
    xp = zp//sigma(zp)
    @assert _padic_index(I, zp * O, p) == (0, 0)
    @assert _padic_index(M, xp * O, p) == (0, 0)
    push!(xps, xp)
  end
  return M, xps, crit_primes
end

function _fractional_ideal_from_base_ring_generators(OE, v)
  # v are elements of the field
  return fractional_ideal(OE, basis_matrix(v) * basis_mat_inv(OE))
end

function _intersect(I::NfRelOrdFracIdl, J::NfRelOrdFracIdl)
  pm = _intersect_modules(basis_pmatrix(I), basis_pmatrix(J))
  return fractional_ideal(order(I), pm)
end

function _intersect_lattice_down_contained(xps, _ps, Lambda)
  _ps_orig = deepcopy(_ps)
  xps_orig = deepcopy(xps)
  if length(_ps) == 0
    return one(nf(Lambda)) * Lambda
  end
  for x in xps
    @assert x in Lambda
  end

  LL = nf(Lambda)(1) * Lambda

  for i in 1:length(xps)
    I = xps[i] * Lambda
    n = norm(I)
    for (p, e) in factor(n)
      if p == _ps[i]
        continue
      end
      if e > 0
        I = inv(p)^(e) * I
      end
    end
    LL = _intersect(LL, I)
  end

  for i in 1:length(xps)
    @assert _padic_index(xps[i] * Lambda, LL, _ps[i]) == (0, 0)
    @assert _padic_index(LL, xps[i] * Lambda, _ps[i]) == (0, 0)
  end

  for q in support(norm(LL))
    if !(q in _ps)
      @assert _padic_index(LL, d * Lambda, q) == (0, 0)
    end
  end

  return LL
end

# Given (x_p)_p, find a fractional ideal I such that I_p = x_p Lambda_p
# for all p and I_q = Lambda_q for all other prime ideals.
#
# The idea is to change Jp = x_p such that (Jp)_p \subseteq Lambda_p
# and Lambda_q \subseteq (Jp)_q for q != p. Then we can just intersect
# the Jp with Lambda.
#
# We first reduce to the case x_p in Lambda for all p by multiplying
# with a common denominator.
function _intersect_lattice_down(xps, _ps, Lambda)
  _ps_orig = copy(_ps)
  xps_orig = copy(xps)

  if length(_ps) == 0
    return one(nf(Lambda)) * Lambda
  end
  d = lcm([denominator(x, Lambda) for x in xps])

  for x in xps
    @assert d*x in Lambda
  end

  _xps = [d * x for x in xps]
  _ps = copy(_ps)

  for p in support(norm(Lambda(d)) * base_ring(Lambda))
    if !(p in _ps)
      push!(_ps, p)
      push!(_xps, nf(Lambda)(d))
    end
  end
  LL = _intersect_lattice_down_contained(_xps, _ps, Lambda)
  pm = deepcopy(basis_pmatrix(LL))
  pm.matrix = divexact(pm.matrix, base_field(nf(Lambda))(d))
  LLL = fractional_ideal(Lambda, pm)

  for q in support(norm(LLL))
    if !(q in _ps)
    end
  end

  for i in 1:length(xps)
    @assert _padic_index(xps_orig[i] * Lambda, LLL, _ps_orig[i]) == (0, 0)
    @assert _padic_index(LLL, xps_orig[i] * Lambda, _ps_orig[i]) == (0, 0)
  end

  for q in support(norm(numerator(LLL)))
    if !(q in _ps_orig)
      @assert _padic_index(LLL, Lambda, q) == (0, 0)
    end
  end

  for q in support(denominator(LLL) * base_ring(Lambda))
    if !(q in _ps_orig)
      @assert _padic_index(LLL, Lambda, q) == (0, 0)
    end
  end

  return LLL
end

function _fractional_ideal_from_base_ring_generators(OE, v, id)
  # v are elements of the field
  return fractional_ideal(OE, pseudo_matrix(basis_matrix(v) * basis_mat_inv(OE), id))
end

function absolute_norm(A::Hecke.AlgAssRelOrdIdl)
  return absolute_norm(norm(A))
end

function absolute_norm(A::Hecke.AlgAssAbsOrdIdl)
  return norm(A)
end

function absolute_norm(A::NfAbsOrdFracIdl)
  return norm(A)
end

function *(a::NfAbsOrdFracIdl{AnticNumberField,nf_elem}, b::AlgAssRelOrdIdl{nf_elem,Hecke.NfAbsOrdFracIdl{AnticNumberField,nf_elem},AlgAss{nf_elem}})
  pm = basis_pmatrix(b)
  pmnew = pseudo_matrix(matrix(pm), map(z -> a * z, coefficient_ideals(pm)))
  return ideal(algebra(order(b)), pmnew)
end

function _padic_index(N, M, p)
  Npb = pseudo_basis(N)
  Mpb = pseudo_basis(M)
  BN = basis_matrix(N)
  BM = basis_matrix(M)
  T = BN * inv(BM)
  n = nrows(BN)
  mini = 0
  for i in 1:n
    for j in 1:n
      if !(T[i, j] == 0)
           mini = min(mini, valuation(T[i, j], p) + valuation(Npb[i][2], p) - valuation(Mpb[i][2], p))
      end
    end
  end
  TT = deepcopy(T)
  pi = elem_in_nf(uniformizer(p))
  for i in 1:n
    for j in 1:n
      TT[i, j] = TT[i, j] * pi^(valuation(Npb[i][2], p) - valuation(Mpb[i][2], p))
    end
  end

  dM = det(basis_pmatrix(M))
  dN = det(basis_pmatrix(N))
  d = dN * inv(dM)
  @assert valuation(d, p) == valuation(det(TT), p)
  return mini, valuation(d, p)
end

function _ambigous_ideals_quotient_basis(OFabs, F, KtoFabs, FabstoF)
  # Given a maximal order OF of an absolute field and K -> F such that F/K is normal, return
  # a basis for {ambiguous ideals of OF with respect to Gal(F/K)}/{ideals of OK}
  res = ideal_type(OFabs)[]
  OF = maximal_order(F)
  dis = discriminant(OF)
  for (p, _) in factor(dis)
    lp = prime_decomposition(maximal_order(F), p)
    c = length(lp)
    @assert length(lp) == 1 && lp[1][2] == 2
    #@show length(c)
    lpabs = prime_decomposition(OFabs, minimum(p))
    phip = 1 * OFabs
    phipp = lp[1][1]
    cc = 0
    for (P, _) in lpabs
      pi = FabstoF(elem_in_nf(p_uniformizer(P)))
      for (Q, _) in lp
        v = valuation(pi, Q)
        @assert 0 <= v <= 1
        if v == 1
          cc += 1
          phip = phip * P
          break
        end
      end
    end
    @assert c == cc
    @assert norm(phip) == absolute_norm(phipp)
    push!(res, phip)
  end
  #@show length(res)
  return res
end

################################################################################
#
#  Binary via ternary
#
################################################################################

function _embed_into_ternary_lattice(L)
  @req ispositive_definite(L) "Lattice must be positive definite"
  V = rational_span(L)
  ML = basis_matrix_of_rational_span(L)
  # To go from V to ambient_space(L), multiply from the right with ML
  # The other direction in general requires solving x * ML = v, if v is an
  # element of rational_span(L)

  G = gram_matrix(V)
  K = base_ring(V)
  Gext = diagonal_matrix(G, matrix(K, 1, 1, [one(K)]))
  W = quadratic_space(K, Gext)
  Lpmat = identity_matrix(K, rank(W))
  OK = base_ring(L)
  coeffs = copy(coefficient_ideals(L))
  push!(coeffs, K(1) * OK)
  Lp = lattice(W, pseudo_matrix(Lpmat, coeffs) )
  return W, Lp, ML
end

function _minimum_of_restriction(L)
  return minimum(_restrict_scalars_absolute(L)[1])
end

function _restriction_represents_one(L)
  return isone(minimum(_restrict_scalars_absolute(L)[1]))
end

function _restrict_scalars_absolute(L)
  V = ambient_space(L)
  E = base_ring(V)
  Vabs, f, g = restrict_scalars(V)
  Babs = absolute_basis(L)
  Babsmat = matrix(E, Babs)
  Mabs = zero_matrix(FlintQQ, length(Babs), rank(Vabs))
  for i in 1:length(Babs)
    v = g(Babs[i])
    for j in 1:length(v)
      Mabs[i, j] = v[j]
    end
  end
  return ZLat(Vabs, Mabs), f, g
end

function _shortest_vectors_restricted(L, m)
  Lre, f, g = _restrict_scalars_absolute(L)
  vecs = short_vectors(Lre, m, m)
  res = []
  for v in vecs
    w = matrix(FlintQQ, 1, rank(Lre), v[1]) * basis_matrix_of_rational_span(Lre)
    push!(res, f(w))
  end
  return res
end

function _orthogonal_complement(M::AbsLat, L::AbsLat)
  V = ambient_space(L)
  M = basis_matrix_of_rational_span(M)
  Morth = orthogonal_complement(V, M)
  # Now intersect KM with L
  pm = _intersection_modules(pseudo_matrix(Morth), pseudo_matrix(L))
  return lattice(V, pm)
end

function _orthogonal_complement(v, L)
  V = ambient_space(L)
  M = matrix(base_ring(V), 1, length(v), v)
  ge = generators(L)
  ge_or = copy(ge)
    for i in 1:length(ge)
    # <v, v> = 1
    ge_or[i] = ge[i] - inner_product(V, ge[i], v) .* v
    @assert inner_product(V, ge_or[i], v) == 0
  end
  pm = pseudo_hnf_kb(pseudo_matrix(transpose(matrix(ge_or))), :lowerleft)
  i = 1
  while iszero_row(pm.matrix, i)
    i += 1
  end

  pm = sub(pm, i:nrows(pm), 1:ncols(pm))

  return lattice(V, pm)
end

function _binary_genus_representatives_via_ternary(L::AbsLat)
  W, Lp, ML = _embed_into_ternary_lattice(L)
  GL = genus(L)
  G = genus_representatives(Lp)
  res = typeof(L)[]
  for LL in G
    vecs = _shortest_vectors_restricted(LL, absolute_degree(nf(base_ring(L))))
    for v in vecs
      LLcomp = _orthogonal_complement(v, LL)
      # Transport this into a full rank lattice
      VV = rational_span(LLcomp)
      LLcomp = lattice(VV, pseudo_matrix(identity_matrix(base_ring(VV), rank(VV)), coefficient_ideals(LLcomp)))
      if genus(LLcomp) == GL
        if any(x -> isisometric(LLcomp, x)[1], res)
          continue
        else
          push!(res, LLcomp)
        end
      end
    end
  end
  return res
end

################################################################################
#
#  Indefinite binary case
#
################################################################################

function _genus_representatives_binary_quadratic_indefinite(_L, max = max, use_auto = true, use_mass = true)
  @assert !isdefinite(_L)
  @assert rank(_L) == 2

  V = ambient_space(_L)
  # 0. Scale to have 1 in Q(V)

  D = diagonal(V)
  _, i = findmin(abs.(norm.(D)))
  d = D[i]

  # so G -> G/d

  L = lattice(quadratic_space(base_ring(V), 1//d * gram_matrix(ambient_space(_L))), pseudo_matrix(_L))

  V = rational_span(L)

  # 1. Find the isometry with (K + K, Tr/2)
  @vprint :GenRep 1 "Determining isometry with CM field ... \n"
  K = base_ring(V)

  mult_tb = Array{elem_type(K), 3}(undef, 2, 2, 2)
  for i in 1:2
    for j in 1:2
      if i != j
        mult_tb[i, j, :] = elem_type(K)[zero(K), zero(K)]
      else
        z = elem_type(K)[zero(K), zero(K)]
        z[i] = one(K)
        mult_tb[i, j, :] = z
      end
    end
  end

  A = AlgAss(K, mult_tb)
  B = basis(A)
  sigma(a) = A([a.coeffs[2], a.coeffs[1]])
  inv2 = inv(A(2))
  phi(x, y) = (x * sigma(y) + y * sigma(x)) * inv2
  G = matrix(K, 2, 2, [phi(B[1], B[1]), phi(B[1], B[2]), phi(B[2], B[1]), phi(B[2], B[2])])
  W = quadratic_space(K, G)
  fl, T = isisometric_with_isometry(V, W)
  @assert fl

  Tinv = inv(T)

  # 2. Transport L to lattice of A.

  pb = pseudo_basis(L)
  # KL has basis a[1] for a in pb
  # So we only need to map the Z-basis of the coefficient ideals
  image_of_first = (T[1, 1] * B[1] + T[1, 2] * B[2])
  image_of_second = (T[2, 1] * B[1] + T[2, 2] * B[2])

  basisofLinA = elem_type(A)[]

  B = absolute_basis(pb[1][2])
  for b in B
    push!(basisofLinA, b * image_of_first)
  end

  B = absolute_basis(pb[2][2])
  for b in B
    push!(basisofLinA, b * image_of_second)
  end

  LinA = ideal_from_lattice_gens(A, basisofLinA)

  Lambda = right_order(LinA)

  OK = maximal_order(K)
  C, mC = class_group(OK)
  cur_class = elem_type(C)[]
  ideal_repr = ideal_type(OK)[]
  p_and_xps = Vector{Tuple{ideal_type(OK), elem_type(A)}}[]
  # We try to find prime ideal representatives
  push!(cur_class, id(C))
  push!(ideal_repr, 1 * OK)
  push!(p_and_xps, Tuple{ideal_type(OK), elem_type(A)}[])
  for p in PrimeIdealsSet(OK, 1, -1)
    if length(cur_class) == order(C)
      break
    end
    c = mC\p
    if !(c in cur_class)
      push!(cur_class, c)
      push!(ideal_repr, p)
      pi = elem_in_nf(uniformizer(p))
      push!(p_and_xps, [(p, A([pi, inv(pi)]))])
    end
  end

  # I have p_and_xps
  _intersect_lattice_down([ p[2] for p in p_and_xps[1]], [ p[1] for p in p_and_xps[1]] , Lambda)

  return Lambda, LinA
end

################################################################################
#
#  Indefinite binary quadratic
#
################################################################################

function _genus_representatives_binary_quadratic_indefinite(L::QuadLat)
  @req degree(base_ring(L)) == 1 "Number field must be of degree 1"
  @req rank(L) == 2 "Lattice must be of rank 2"
  K = nf(base_ring(L))
  f, e = _lattice_to_binary_quadratic_form(L) # e is the scaling factor
  d = discriminant(f)
  @assert d > 0
  cls = _equivalence_classes_binary_quadratic_indefinite(d, proper = false, primitive = false)
  res = typeof(f)[f]
  K, = rationals_as_number_field()
  G = genus(_binary_quadratic_form_to_lattice(f, K, e))
  lat = typeof(L)[]

  for g in cls
    LL = _binary_quadratic_form_to_lattice(g, K, e)
    GG = genus(LL)
    if GG == G
      push!(lat, LL)
    end
  end
  return lat
end

function _lattice_to_binary_quadratic_form(L::QuadLat)
  M = absolute_basis_matrix(L) # This corresponds to a basis of L
  @assert nrows(M) == 2 && ncols(M) == 2
  G = gram_matrix(ambient_space(L), M)
  GG = change_base_ring(FlintQQ, G)
  d = denominator(GG)
  f = binary_quadratic_form(FlintZZ(d * GG[1, 1]), FlintZZ(2 * d * GG[1, 2]), FlintZZ(d * GG[2, 2]))
  return f, d
end

function _equivalence_classes_binary_quadratic_indefinite(d::fmpz; proper::Bool = false, primitive::Bool = true)
  if primitive
    return _equivalence_classes_binary_quadratic_indefinite_primitive(d, proper = proper)
  else
    res = QuadBin{fmpz}[]
    for n in Divisors(d, units = false, power = 2) # n^2 | d
      cls = _equivalence_classes_binary_quadratic_indefinite_primitive(divexact(d, n^2), proper = proper)
      for f in cls
        push!(res, n*f)
      end
    end
    return res
  end
end

function _equivalence_classes_binary_quadratic_indefinite_primitive(d::fmpz; proper::Bool = false)
  @assert d > 0
  Qx = Hecke.Globals.Qx
  x = gen(Qx)
  f = x^2 - d * x + (d^2 - d)//4
  @assert isone(denominator(f))
  K, a = number_field(f, "a", cached = false) # a is (d + \sqrt(d))//2
  O = equation_order(K)
  C, _dlog, _exp = narrow_picard_group(O)
  res = QuadBin{fmpz}[]
  # C gives me all proper classes of definit forms
  # So if proper = true, we don't have to do anything
  # and if proper = false, we have to sieve using isequivalent
  for c in C
    I = _exp(c)
    J = numerator(I)
    f = _ideal_to_form(J, d)
    if proper || all(h -> !isequivalent(h, f, proper = false), res)
      push!(res, reduction(f))
    end
  end
  return res
end

# This is from Kani
function _ideal_to_form(I::NfAbsOrdIdl, delta)
  # first make primitive
  M = _hnf(basis_matrix(I), :lowerleft)
  g = reduce(gcd, [M[1, 1], M[1, 2], M[2, 2]])
  M = divexact(M, g)
  B = M[2, 1]
  C = M[2, 2]
  a = M[1, 1]
  @assert isone(C)
  b = -(2 * B + delta)
  c = -divexact(divexact(delta - b^2, 4), a)
  @assert -4 * a * c == delta - b^2
  f = binary_quadratic_form(a, b, c)
  @assert discriminant(f) == delta
  return f
end

function automorphism_group_generators(g::QuadBin{fmpz})
  d = discriminant(g)
  @assert d > 0
  Qx = Hecke.Globals.Qx
  x = gen(Qx)
  f = x^2 - d * x + (d^2 - d)//4
  @assert isone(denominator(f))
  K, a = number_field(f, "a", cached = false) # a is (d + \sqrt(d))//2
  O = equation_order(K)
  deltasqrt = O(2 * a - discriminant(f))
  @assert deltasqrt^2 == d

  U, mU = unit_group(O)
  A = abelian_group([2])
  h = hom(U, A, [norm(mU(U[1])) == -1 ? A([1]) : A([0]), norm(mU(U[2])) == -1 ? A([1]) : A([0])])
  k, mk = kernel(h, false)
  @assert ngens(k) == 2
  norm_one_gens = elem_type(O)[mU((mk(k[1]))), mU((mk(k[2])))]
  # We need to write the elements in the basis 1//2, deltasqrt//2
  # O has basis 1, (d + sqrt(d)//2)
  gens = dense_matrix_type(FlintZZ)[]
  for i in 1:2
    a, b = coordinates(norm_one_gens[i])
    _x, _y = (2*a + b * d, b)
    @assert norm_one_gens[i] == divexact(_x + _y * deltasqrt, 2)
    @assert iseven(_x + g.b * _y)
    @assert iseven(_x + d * _y)
    T = zero_matrix(FlintZZ, 2, 2)
    T[1, 1] = divexact(_x + _y * g.b, 2)
    T[1, 2] = g.c * _y
    T[2, 1] = -g.a * _y
    T[2, 2] = divexact(_x - _y * g.b, 2)
    push!(gens, T)
  end
  # Now test if g is ambiguous or not
  gg = binary_quadratic_form(g.a, -g.b, g.c)
  fl = isequivalent(g, gg, proper = true)
  if fl
    fll, q = haspreimage(h, A([1]))
    elem = mU(q)
    a, b = coordinates(elem)
    _x, _y = (2*a + b * d, b)
    @assert elem == divexact(_x + _y * deltasqrt, 2)
    @assert iseven(_x + g.b * _y)
    @assert iseven(_x + d * _y)
    T = zero_matrix(FlintZZ, 2, 2)
    T[1, 1] = divexact(_x + _y * g.b, 2)
    T[1, 2] = g.c * _y
    T[2, 1] = -g.a * _y
    T[2, 2] = divexact(_x - _y * g.b, 2)
    TT = matrix(FlintZZ, 2, 2, [0, -1, 1, 0])
    @assert det(TT * T) == -1
    push!(gens, TT * T)
  end
  @assert all(T -> _action(g, T) == g, gens)
  return gens
end

function _form_to_ideal(f::QuadBin{fmpz}, O, a)
  # a must be d + sqrt(d)//2 and O = ZZ[a]
  deltasqrt = O(2 * a - discriminant(f))
  # deltasqrt^2 == delta
  @assert deltasqrt^2 == discriminant(f)
  _a = f[1]
  _b = f[2]
  _c = f[3]
  return ideal(O, [O(_a), divexact(O(-_b + deltasqrt), 2)])
end

function _genus_representatives_binary_quadratic_indefinite(f::QuadBin{fmpz})
  d = discriminant(f)
  @assert d > 0
  cls = _equivalence_classes_binary_quadratic_indefinite(d, proper = false)
  res = typeof(f)[]
  K, = rationals_as_number_field()
  L = _binary_quadratic_form_to_lattice(f, K)
  G = genus(L)
  for g in cls
    GG = genus(_binary_quadratic_form_to_lattice(f, K))
    if G == GG
      push!(res, g)
    end
  end
  return res
end

function _binary_quadratic_form_to_lattice(f::QuadBin{fmpz}, K, e::fmpz = fmpz(1))
  a = f[1]
  b = f[2]
  c = f[3]
  G = matrix(K, 2, 2, [a//(e), b//(2*e), b//(2*e), c//e])
  L = lattice(quadratic_space(K, G), identity_matrix(K, 2))
end

function _equivalence_classes_binary_quadratic_form_reducible(d::fmpz; proper::Bool = false, primitive::Bool = true)
  if primitive
    return _equivalence_classes_binary_quadratic_form_reducible_primitive(d, proper = proper)
  else
    res = QuadBin{fmpz}[]
    for n in Divisors(d, units = false, power = 2) # n^2 | d
      cls = _equivalence_classes_binary_quadratic_form_reducible_primitive(divexact(d, n^2), proper = proper)
      for f in cls
        push!(res, n*f)
      end
    end
    return res
  end

end

function _equivalence_classes_binary_quadratic_form_reducible_primitive(N::fmpz; proper::Bool = false)
  d = sqrt(N)
  @assert d^2 == N
  res = fmpz[]
  for i in 1:(d-1)
    if !isone(gcd(i, d))
      continue
    end
    if proper
      push!(res, i)
    else
      j = invmod(i, d)
      if j in res
        continue
      end
      push!(res, i)
    end
  end
  return QuadBin{fmpz}[binary_quadratic_form(i, d, zero(fmpz)) for i in res]
end


text:
using MatchingMarkets.Util

@testset "Testing util.jl" begin

    @testset "binary max heap" begin
        bh1 = Util.BinMaxHeap(10)
        bh2 = Util.BinMaxHeap(10)
        bh3 = Util.BinMaxHeap(10)
        bh4 = Util.BinMaxHeap(10)
        bh5 = Util.BinMaxHeap(10)

        @testset "empty heap equality" begin
            @test bh1 == bh2
        end

        @testset "accessing an empty heap" begin
            @test_throws BoundsError Util.top(bh1)
        end

        @testset "pop! from an empty heap" begin
            @test_throws BoundsError Util.pop!(bh1)
        end

        @testset "heap ordered" begin
            for i in [3, 2, 10, 6, 5, 9, 7, 1, 4, 8]
                Util.push!(bh2, i)
            end
            for i in 10:-1:1
                @test Util.pop!(bh2) == i
            end
            for i in [1, 8, 2, 3, 6, 5, 4, 7, 10, 9]
                Util.push!(bh2, i)
            end
            for i in 10:-1:1
                @test Util.pop!(bh2) == i
            end
        end

        @testset "length of heap" begin
            for i in 1:10
                Util.push!(bh3, i)
                @test length(bh3) == i
            end
        end

        @testset "pushing to full heap not allowed" begin
            for i in 1:10
                Util.push!(bh4, i)
            end
            @test_throws BoundsError Util.push!(bh4, 11)
        end

        @testset "least replacement" begin
            for i in 10:-1:1
                Util.push!(bh5, i)
            end
            Util.replace_least!(bh5, 11)
            @test Util.top(bh5) == 11
            Util.replace_least!(bh5, 10)
            @test Util.top(bh5) == 10
            Util.replace_least!(bh5, 9)
            @test Util.top(bh5) == 9
        end
    end

end

ext: jl


text:
function scale_polytope!(model::MetNet, factor::Real)

    for (rxni, rxn) in enumerate(model.rxns)
        lb, ub = bounds(model, rxni)
        lb -= abs(lb) * factor
        ub += abs(ub) * factor
        bounds!(model, rxni, lb, ub)
    end
    return model
end


text:
using MOCNeutronTransport

# I can't remember if the bottom left corner of the CMFD grid needs to be (0, 0, 0)
# To compensate, we do some extra math to ensure that happens.
grid_offset = 0.1 # Buffer between boundary of the model and the edge of the CMFD grid (cm)
# Bounding box for the CMFD grid
# xmin, xmax, ymin, ymax
bb = (0.0, 7.46252 + 2*grid_offset,
      0.0, 8.31535 + 2*grid_offset)

# Model setup
# ----------------------------------------------------------------------------------------
gmsh.initialize()
# Import CAD file
gmsh.merge("xsec.step")
# Visualize the model
# gmsh.fltk.run()

# Shift the model to the origin of the xy-plane
# Get the dim tags of all dimension 2 entites
dim_tags = gmsh.model.get_entities(2)
# Rotate the entities in dim_tags in the built-in CAD representation by angle radians
# around the axis of revolution defined by the point (x, y, z) and the direction (ax, ay, az).
# gmsh.model.occ.rotate(dim_tags, x, y, z, ax, ay, az, angle)
gmsh.model.occ.rotate(dim_tags, 0, 1, 0, 0, 1, 0, -π/2)
# Translate it so that the corner is at the origin
gmsh.model.occ.translate(dim_tags, 37.3225 - 0.0099 + 10*grid_offset,
                                   44.09 - 0.00174209 + 10*grid_offset,
                                   -436.5625)
# Model is in mm. Convert to cm by shrinking everything by 1/10
gmsh.model.occ.dilate(dim_tags, 0, 0, 0, 1//10, 1//10, 0)
# Synchronize the CAD and Gmsh models
gmsh.model.occ.synchronize()
# Verify results
# gmsh.fltk.run()

# Assign materials
# ------------------------------------------------------------------------------------------
# Elements containing fuel were determined visually, since the numbering scheme isn't obvious
fuel_tags = [4, 5, 6, 9, 10, 11, 12, 13, 15, 20, 21, 25, 27, 31, 32, 33, 36, 38]
clad_tags = [t[2] for t in dim_tags]
# Determine which tags are clad by doing fuel\clad
filter!(x->x ∉  fuel_tags, clad_tags)
# Assign the materials as physical groups
p = gmsh.model.add_physical_group(2, fuel_tags)
gmsh.model.set_physical_name(2, p, "MATERIAL_UO2")
p = gmsh.model.add_physical_group(2, clad_tags)
gmsh.model.set_physical_name(2, p, "MATERIAL_CLAD")
# Verify results
# gmsh.fltk.run()

# Overlay CMFD/hierarchical grid
# -------------------------------------------------------------------------------------------
# The arguments to gmsh_overlay_rectangular_grid are:
# (bb, material, nx, ny) or (bb, material, x, y)
#
# bb: The bounding box for the grid
#   bb = (xmin, xmax, ymin, ymax)::NTuple{4, AbstractFloat}
# material: The material with which to fill empty space/give to entities without a material
#   material::String
#   Note the material string must start with something that when converted to uppercase
#       yields "MATERIAL_", ex: "Material_A" or "material_A"
# (nx, ny)/(x, y): The hierarchical grid divisions.
#   The grid is mapped to MPACT in the following manner:
#       1st level: Assembly map
#       2nd level: Assemblies/lattices (same map, since 2D)
#       3rd level: Ray tracing modules (must be equal size)
#       4th level: Coarse cells in the CMFD grid
#   This has two input variations
#   Input 1) nx, ny
#       The number of equal x and y divisions in each level of the hierarchy.
#       nx::Vector{Int}, ny::Vector{Int}
#       Example:
#           nx = [3, 2, 1, 1]
#           The 1st level (assembly map) will be split into 3 pieces in x.
#           The 2nd level will split each of the 1st level grid elements into 2 piece in x
#               Theses elements are now 1/6 the total grid's x width
#           The 3rd level will be the same as above
#           The 4th level will be the same as above
#   Input 2) x, y
#       The x and y locations of the divisions. These need not be equally spaced.
#       x::Vector{Vector{T}}, y::Vector{Vector{T}}
#       Note that end points don't need to be included, and that any divisions in the coarser
#           grid levels are propagated downward.
#       Example:
#          For bb = (0, 1, 0, 1)
#          x = [[0.5],[0.2, 0.75]]
#          is equal to
#          x = [[0.0, 0.5, 1.0], [0.0, 0.2, 0.5, 0.75, 1.0]]
#   Final tidbit:
#       Any missing enities are assumed to be 1 division, where the last element in the input
#           vector is assumed to be the coarse cells.
#       Example:
#           nx = [1, 1, 4, 3] gives the same output as [4, 3]
#           4 RT modules (in x), 3 coarse cells (in x) per module = 12 total coarse cells (in x)
nx = [9, 1]
ny = [10, 1]
grid_tags = gmsh_overlay_rectangular_grid(bb, "MATERIAL_MODERATOR", nx, ny)
# Verify results
# gmsh.fltk.run()

# Mesh
# ------------------------------------------------------------------------------------------------
# Set the characteristic edge length of the mesh cells
lc = 0.3 # cm
gmsh.model.mesh.set_size(gmsh.model.get_entities(0), lc)
# Optional mesh optimization:
niter = 2 # The optimization iterations

# Triangles
# gmsh.model.mesh.generate(2) # 2 is dimension of mesh
# for () in 1:niter
#     gmsh.model.mesh.optimize("Laplace2D")
# end

# Quadrilaterals
# The default recombination algorithm might leave some triangles in the mesh, if
# recombining all the triangles leads to badly shaped quads. In such cases, to
# generate full-quad meshes, you can either subdivide the resulting hybrid mesh
# (with `Mesh.SubdivisionAlgorithm' set to 1), or use the full-quad
# recombination algorithm, which will automatically perform a coarser mesh
# followed by recombination, smoothing and subdivision.
# Mesh recombination algorithm (0: simple, 1: blossom, 2: simple full-quad, 3: blossom full-quad)
# Default = 1
# gmsh.option.set_number("Mesh.RecombineAll", 1) # recombine all triangles
# gmsh.option.set_number("Mesh.Algorithm", 8) # Frontal-Delaunay for quads. Better 2D algorithm
# gmsh.option.set_number("Mesh.RecombinationAlgorithm", 1)
# gmsh.model.mesh.generate(2)
# for () in 1:niter
#     gmsh.model.mesh.optimize("Laplace2D")
#     gmsh.model.mesh.optimize("Relocate2D")
#     gmsh.model.mesh.optimize("Laplace2D")
# end

# 2nd order triangles
gmsh.model.mesh.generate(2) # Triangles first for high order meshes.
gmsh.model.mesh.set_order(2)
for () in 1:niter
    gmsh.model.mesh.optimize("HighOrderElastic")
    gmsh.model.mesh.optimize("Relocate2D")
    gmsh.model.mesh.optimize("HighOrderElastic")
end

# 2nd order quadrilaterals
# These can be problematic for large lc. They have trouble respecting CAD boundaries.
# gmsh.option.set_number("Mesh.RecombineAll", 1) # recombine all triangles
# gmsh.option.set_number("Mesh.Algorithm", 8) # Frontal-Delaunay for quads. Better 2D algorithm
# gmsh.option.set_number("Mesh.RecombinationAlgorithm", 3)
# gmsh.model.mesh.generate(2)
# gmsh.model.mesh.set_order(2)
# for () in 1:niter
#     gmsh.model.mesh.optimize("HighOrderElastic")
#     gmsh.model.mesh.optimize("Relocate2D")
#     gmsh.model.mesh.optimize("HighOrderElastic")
# end

# gmsh.fltk.run()

# Mesh conversion
# ----------------------------------------------------------------------------------------------------
# We want to write an Abaqus file from Gmsh, read the Abaqus file into Julia, and convert that into
# a rectangularly partitioned hierarchical XDMF file, the final input to MPACT
#
# Write Abaqus file
gmsh.write("fnr.inp")
gmsh.finalize() # done with Gmsh. Finalize
mesh = read_abaqus_2d("fnr.inp")
# write_vtk_2d("fnr.vtk", mesh) # Can print vtk for sanity check. No cell sets
# write_xdmf_2d("fnr.xdmf", mesh) # Can print single grid XDMF
# Convert mesh into Hierarchical Rectangularly Partitioned Mesh
HRPM = partition_rectangularly(mesh)
write_xdmf_2d("fnr.xdmf", HRPM)

# Mass conservation
# ----------------------------------------------------------------------------------------------------
# See how well your mesh conserves mass. Analytic solution not available, so fine tri6 mesh used as
# reference
fuel_area = area(mesh, "MATERIAL_UO2")
fuel_area_ref = 6.119632103711788
println("Fuel area (reference): $fuel_area_ref")
println("Fuel area ( computed): $fuel_area")
err = 100*(fuel_area - fuel_area_ref)/fuel_area_ref
println("Error (relative): $err %")
println("")

clad_area = area(mesh, "MATERIAL_CLAD")
clad_area_ref = 19.884367501087866
println("Clad area (reference): $clad_area_ref")
println("Clad area ( computed): $clad_area")
err = 100*(clad_area - clad_area_ref)/clad_area_ref
println("Error (relative): $err %")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using IntervalArithmetic
if VERSION < v"0.7.0-DEV.2004"
    using Base.Test
else
    using Test
end

@testset "Interval loop tests" begin
    i = 1

    @test Interval(i,i).lo == 1
    @test @interval(i).lo == 1

    for i in 1:10
        a = @interval(i)
        @test a.lo == i
    end

end

## Calculate pi by summing 1/i^2 to give pi^2/6:

setprecision(Interval, 53)

function calc_pi1(N)
    S1 = @interval(0)

    for i in 1:N
        S1 += @interval(1/i^2)
    end
    S1 += @interval(1/(N+1), 1/N)

    sqrt(6*S1)
end

function calc_pi2(N)
    S2 = @interval(0)

    for i in 1:N
        S2 += 1/i^2
    end
    S2 += @interval(1/(N+1), 1/N)

    sqrt(6*S2)
end

function calc_pi3(N)
    S3 = @floatinterval(0)

    for i in 1:N
        S3 += 1/i^2
    end
    S3 += @floatinterval(1/(N+1), 1/N)

    sqrt(6*S3)
end

function calc_pi4(N)
    S4 = @floatinterval(0)
    II = @floatinterval(1)

    for i in N:-1:1
        S4 += II / (i^2)
    end
    S4 += II / @floatinterval(N, N+1)

    sqrt(6*S4)
end

function calc_pi5(N)
    S5 = @floatinterval(0)
    II = @floatinterval(1)

    for i in N:-1:1
        S5 += 1 // (i^2)
    end
    S5 += 1 / @floatinterval(N, N+1)

    sqrt(6*S5)
end

@testset "Pi tests" begin

    big_pi = setprecision(256) do
        big(pi)
    end

    N = 10000
    pi1 = calc_pi1(N)
    pi2 = calc_pi2(N)
    pi3 = calc_pi3(N)
    pi4 = calc_pi4(N)
    pi5 = calc_pi5(N)

    @test big_pi ∈ pi1
    @test big_pi ∈ pi2
    @test big_pi ∈ pi3
    @test big_pi ∈ pi4
    @test big_pi ∈ pi5

    @test pi1 == pi2
    @test pi2 == pi3

end

setprecision(Interval, Float64)


text:
module Serialization
    using SparseArrays

    export write_sparse_matrix
    export read_sparse_matrix

    # Writes the sparse matrix M in text format to a file with the given name
    function write_sparse_matrix(filename, M, T=Int64)
        open(filename, "w") do f
            write(f, "$(T(size(M, 1)))\n$(T(size(M, 2)))\n")
            nz = findnz(M)
            for n = 1:2
                foreach(x -> write(f, "$(T(x))\n"), nz[n])
            end
            foreach(x -> write(f, "$(Float64(x))\n"), nz[3])
        end
    end

    # Reads and returns a sparse matrix written by write_sparse_matrix
    function read_sparse_matrix(filename, T=Int64)
        data = open(filename) do f
            readlines(f)
        end
        m = parse(T, data[1])
        n = parse(T, data[2])
        c = Integer((length(data) - 2) / 3)
        I = map(x -> parse(T, x), data[3:3+c-1])
        J = map(x -> parse(T, x), data[3+c:3+2*c-1])
        V = map(x -> parse(Float64, x), data[3+2*c:3+3*c-1])
        return sparse(I, J, V, m, n)
    end
end


text:
using LaTeXStrings, Compat
using Compat.Test

tst1 = L"an equation: $\alpha^2$"
tst1s = String(tst1)
@test tst1 == "an equation: \$\\alpha^2\$" == tst1s
@test L"\alpha^2" == "\$\\alpha^2\$"
@test repr("text/latex", tst1) == tst1s
@test sprint(show, "text/latex", tst1) == tst1s
@test latexstring(tst1s) == tst1 == LaTeXString(tst1s)

@test ccall(:strlen, Csize_t, (Cstring,), tst1) == ccall(:strlen, Csize_t, (Ptr{UInt8},), tst1) == sizeof(tst1)

for f in (length, eachindex, pointer, collect, ncodeunits, codeunits)
    @test f(tst1) == f(tst1.s)
end
@test nextind(tst1, 1) == prevind(tst1, 3) == 2
@test isvalid(tst1, 1)
@test findnext(isequal('o'), L"foo", 1) == 3
@test tst1[1] == 'a'
@test codeunit(tst1, 1) == UInt8('a')
@test L"foo"[1] == '$'
@test tst1[1:2] == tst1[0x01:0x02] == tst1[[1,2]] == "an"


text:
"""
   lq_estimation_data()

Returns a 2 x 100 matrix with the measurement data for the linear quadratic
state estimation example of section 17.3.1.
"""
function lq_estimation_data()
# data for linear quadratic estimation example
return ymeas = [ 1.70778099e+01 -1.83924162e+01 6.10394456e+00 -7.97507214e+00 2.27355721e+00 2.34308942e+01 3.39867360e+01 4.24050167e+01 2.05731680e+01 4.94395071e+01 3.45442498e+01 7.88341002e+01 8.67009078e+01 7.22735074e+01 2.31850194e+01 1.25784108e+02 1.11441506e+02 1.34603580e+02 9.08335516e+01 1.39211826e+02 1.82597858e+02 1.55620959e+02 1.35394535e+02 1.71816643e+02 1.79738746e+02 2.62311327e+02 1.65873945e+02 1.56867751e+02 2.23467712e+02 2.53129175e+02 2.03679219e+02 2.20062198e+02 2.81625929e+02 2.38298096e+02 3.01005533e+02 2.74101690e+02 2.41510019e+02 2.48551509e+02 2.25600934e+02 1.54680964e+02 1.98224063e+02 2.03866603e+02 2.08089219e+02 2.08244771e+02 1.69046465e+02 1.74862089e+02 1.93984500e+02 1.61771170e+02 1.44928387e+02 1.07953687e+02 1.42606362e+02 1.01143225e+02 1.93119878e+02 4.86121107e+01 1.40429785e+02 1.16747984e+02 7.20067571e+01 3.33845460e+01 1.73703911e+00 4.13324794e+01 5.10008943e+01 7.78212231e+01 9.73475231e+01 7.20156638e+00 4.18824911e+01 3.73829417e+01 1.17537300e+01 6.50543261e+01 7.01591558e+00 3.89538168e+01 1.04449893e+02 6.21834962e+01 9.85315111e+01 5.14442316e+01 4.75120937e+01 6.13570702e+01 5.97455584e+01 4.39494903e+01 8.26649843e+01 6.08501125e+01 1.07528400e+02 1.62298903e+02 1.60557305e+02 1.40797792e+02 1.39902849e+02 1.17845958e+02 2.00958386e+02 2.02880814e+02 1.76931913e+02 1.65895483e+02 2.01597849e+02 1.98339869e+02 2.17976883e+02 2.61641268e+02 2.04828587e+02 2.84764612e+02 2.22338022e+02 2.95402558e+02 2.60828319e+02 2.24987351e+02 ;
 1.63811511e+01 -1.02102123e+00 -2.62172723e+00 1.70421017e+01 -8.77394913e+00 -5.20869630e+01 -1.60324186e+01 -4.24246849e+01 -4.91794661e+01 -3.06611440e+01 -1.52727993e+01 -8.11660285e+01 -1.16511615e+02 -1.41130660e+01 -4.61581740e+01 -7.50780171e+01 -1.18182293e+02 -1.29205129e+02 -1.64972334e+02 -1.02577288e+02 -1.22567633e+02 -1.81248284e+02 -1.27278951e+02 -2.04119652e+02 -2.27687535e+02 -2.51393552e+02 -2.74520269e+02 -2.66083141e+02 -3.11282088e+02 -3.05291240e+02 -3.87589339e+02 -4.13509577e+02 -3.79896889e+02 -4.66118373e+02 -4.44740533e+02 -4.57258956e+02 -4.92460236e+02 -4.94789750e+02 -5.42765233e+02 -5.26015331e+02 -5.90580789e+02 -6.21374263e+02 -6.17038815e+02 -6.21125522e+02 -6.39106527e+02 -6.94088028e+02 -7.14041092e+02 -7.35457565e+02 -7.26167922e+02 -7.74166470e+02 -8.35004328e+02 -8.56688388e+02 -8.68452776e+02 -9.07279418e+02 -8.85158459e+02 -9.42174345e+02 -9.96411241e+02 -9.31762579e+02 -1.01901698e+03 -1.09368046e+03 -9.92186943e+02 -1.10887128e+03 -1.11942878e+03 -1.18326208e+03 -1.18145562e+03 -1.14747530e+03 -1.19375816e+03 -1.22016814e+03 -1.22906150e+03 -1.28998644e+03 -1.24880527e+03 -1.30671968e+03 -1.36047927e+03 -1.39128750e+03 -1.37220909e+03 -1.39419928e+03 -1.44482396e+03 -1.47010103e+03 -1.42228170e+03 -1.44700843e+03 -1.49317030e+03 -1.55328302e+03 -1.47098557e+03 -1.54946659e+03 -1.55488173e+03 -1.54173960e+03 -1.56711254e+03 -1.59185254e+03 -1.56509343e+03 -1.60766101e+03 -1.65272075e+03 -1.65834577e+03 -1.56199411e+03 -1.64068671e+03 -1.62112540e+03 -1.63171892e+03 -1.66863768e+03 -1.63170753e+03 -1.65546697e+03 -1.64358222e+03
]

end

max_issues_repo_name: lynnlangit/julia-linear-algebra


text:
using DiffRules, SpecialFunctions, NaNMath
using Base.FastMath: fast_op, make_fastmath

@nograd isinf, isnan, isfinite, div

# TODO use CSE here

for (M, f, arity) in DiffRules.diffrules()
  arity == 1 || continue
  Δ = :Δ
  dx = DiffRules.diffrule(M, f, :x)
  if f in [:abs, :abs2]
    Δ = :(real($Δ))
  else
    dx = :(conj($dx))
  end
  @eval begin
    @adjoint $M.$f(x::Number) = $M.$f(x),
      Δ -> ($Δ * $dx,)
  end
end

for (M, f, arity) in DiffRules.diffrules()
  arity == 2 || continue
  f == :^ && continue
  da, db = DiffRules.diffrule(M, f, :a, :b)
  @eval begin
    @adjoint $M.$f(a::Number, b::Number) = $M.$f(a, b),
      Δ -> (Δ * conj($da), Δ * conj($db))
  end
end

@adjoint Base.:^(x::Number, p::Number) = x^p,
  Δ -> (Δ * conj(p * x^(p-1)), Δ * conj(x^p * log(complex(x))))
@adjoint Base.literal_pow(::typeof(^), x::Number, ::Val{p}) where {p} =
  Base.literal_pow(^,x,Val(p)),
  Δ -> (nothing, Δ * conj(p * Base.literal_pow(^,x,Val(p-1))), nothing)

@adjoint Base.convert(T::Type{<:Real}, x::Real) = convert(T, x), ȳ -> (nothing, ȳ)
@adjoint (T::Type{<:Real})(x::Real) = T(x), ȳ -> (nothing, ȳ)

for T in Base.uniontypes(Core.BuiltinInts)
    @adjoint (::Type{T})(x::Core.BuiltinInts) = T(x), Δ -> (Δ,)
end

@adjoint Base.:+(xs::Number...) = +(xs...), Δ -> map(_ -> Δ, xs)

@adjoint Base.muladd(x::Number, y::Number, z::Number) =
  Base.muladd(x, y, z), ō -> (y'ō, x'ō, ō)

@adjoint Base.fma(x::Number, y::Number, z::Number) =
  Base.fma(x, y, z), ō -> (y'ō, x'ō, ō)

@adjoint function sincos(x)
  s, c = sincos(x)
  (s, c), ((s̄, c̄),) -> (s̄*c - c̄*s,)
end

@adjoint acosh(x::Complex) =
  acosh(x), Δ -> (Δ * conj(inv(sqrt(x - 1) * sqrt(x + 1))),)

@adjoint a // b = (a // b, c̄ -> (c̄ * 1//b, - c̄ * a // b // b))

@nograd floor, ceil, trunc, round, hash

# Complex Numbers

@adjoint (T::Type{<:Complex})(re, im) = T(re, im), c̄ -> (nothing, real(c̄), imag(c̄))

@adjoint real(x::Number) = real(x), r̄ -> (real(r̄),)
@adjoint conj(x::Number) = conj(x), r̄ -> (conj(r̄),)
@adjoint imag(x::Number) = imag(x), ī -> (real(ī)*im,)

DiffRules._abs_deriv(x::Complex) = x/abs(x)

 # adjoint for Fastmath operations
for (f, fastf) in fast_op
  if DiffRules.hasdiffrule(:Base, f, 1)
    dx = DiffRules.diffrule(:Base, f, :x)
    Δ = :Δ
    if f in [:abs, :abs2]
      Δ = :(real($Δ))
    else
      dx = :(conj($dx))
    end
    @eval begin
      @adjoint Base.FastMath.$fastf(x::Number) =
        Base.FastMath.$fastf(x), Δ -> ($Δ * make_fastmath($dx),)
    end
  elseif DiffRules.hasdiffrule(:Base, f, 2)
    dx, dy = DiffRules.diffrule(:Base, f, :x, :y)
    @eval begin
      @adjoint Base.FastMath.$fastf(x::Number, y::Number) =
        Base.FastMath.$fastf(x, y),
        Δ -> (Δ * make_fastmath(conj($dx)), Δ * make_fastmath(conj($dy)))
    end
  end
end


text:
function compute()
N=1000
ans=0
for i=1:N-1
	if(i%3==0|| i%5==0)
		global ans+=i
	end
end
ans
end

print(compute())


text:
# point = Point(locaslID, x, y, left, right, flag_1, flag_2, nbhs, conn, nx, ny, prim, flux_res, q, dq, entropy, xpos_nbhs, xneg_nbhs, ypos_nbhs, yneg_nbhs, xpos_conn, xneg_conn, ypos_conn, yneg_conn, delta)

function setNormals(self::Point, n)
    self.nx = n[1]
    self.ny = n[2]
    return nothing
end

function getxy(self)
    return (self.x, self.y)
end

function setConnectivity(self, conn)
    self.xpos_conn = conn[1]
    self.xpos_nbhs = length(conn[1])
    self.xneg_conn = conn[2]
    self.xneg_nbhs = length(conn[2])
    self.ypos_conn = conn[3]
    self.ypos_nbhs = length(conn[3])
    self.yneg_conn = conn[4]
    self.yneg_nbhs = length(conn[4])
    return nothing
end

# function setSmallestPointDistance(self::Point, distance)
#     self.short_distance = distance
# end

function convertToFixedArray(targetArray1, originalStruct::Point, idx)
    targetArray1[idx] = FixedPoint(originalStruct.localID,
                                        originalStruct.x,
                                        originalStruct.y,
                                        originalStruct.left,
                                        originalStruct.right,
                                        originalStruct.flag_1,
                                        originalStruct.flag_2,
                                        originalStruct.short_distance,
                                        originalStruct.nbhs,
                                        originalStruct.nx,
                                        originalStruct.ny,
                                        0.0
                                    )
    return nothing
end

function convertToNeighbourArray(targetArray2, originalStruct::Point, idx)
    targetArray2[1] = originalStruct.xpos_nbhs
    targetArray2[2] = originalStruct.xneg_nbhs
    targetArray2[3] = originalStruct.ypos_nbhs
    targetArray2[4] = originalStruct.yneg_nbhs
    targetArray2[5:4 + originalStruct.nbhs, idx] = originalStruct.conn
    targetArray2[15:14 + originalStruct.xpos_nbhs, idx] = originalStruct.xpos_conn
    targetArray2[25:24 + originalStruct.xneg_nbhs, idx] = originalStruct.xneg_conn
    targetArray2[35:34 + originalStruct.ypos_nbhs, idx] = originalStruct.ypos_conn
    targetArray2[45:44 + originalStruct.yneg_nbhs, idx] = originalStruct.yneg_conn
    return nothing
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver