text:
export genus, representative, rank, det, uniformizer, det_representative,
       gram_matrix, representative, genus, genera_hermitian,
       local_genera_hermitian, rank, orthogonal_sum

################################################################################
#
#  Local genus symbol
#
################################################################################

# Need to make this type stable once we have settled on a design
mutable struct LocalGenusHerm{S, T}
  E::S                                # Field
  p::T                                # prime of base_field(E)
  data::Vector{Tuple{Int, Int, Int}}  # data
  norm_val::Vector{Int}               # additional norm valuation
                                      # (for the dyadic case)
  isdyadic::Bool                      # 2 in p
  isramified::Bool                    # p ramified in E
  non_norm_rep                        # u in K*\N(E*)
  ni::Vector{Int}                     # ni for the ramified, dyadic case

  function LocalGenusHerm{S, T}() where {S, T}
    z = new{S, T}()
    return z
  end
end

local_genus_herm_type(E) = LocalGenusHerm{typeof(E), ideal_type(order_type(base_field(E)))}

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

function Base.show(io::IO, ::MIME"text/plain", G::LocalGenusHerm)
  compact = get(io, :compact, false)
  if !compact
    if isdyadic(G) && isramified(G)
       print(io, "Local genus symbol (scale, rank, det, norm) at ")
    else
      print(io, "Local genus symbol (scale, rank, det) at ")
    end
    if length(G) > 0
      print(IOContext(io, :compact => true), prime(G), ":")
      print(io, "\n")
    else
      print(IOContext(io, :compact => true), prime(G), " of rank zero")
    end
  end
  if isdyadic(G) && isramified(G)
    for i in 1:length(G)
      print(io, "(", scale(G, i), ", ", rank(G, i), ", ",
            det(G, i) == 1 ? "+" : "-", ", ", norm(G, i), ")")
    end
  else
    for i in 1:length(G)
      print(io, "(", scale(G, i), ", ", rank(G, i), ", ",
            det(G, i) == 1 ? "+" : "-",  ")")
    end
  end
end

function Base.show(io::IO, G::LocalGenusHerm)
  if isdyadic(G) && isramified(G)
    for i in 1:length(G)
      print(io, "(", scale(G, i), ", ", rank(G, i), ", ",
            det(G, i) == 1 ? "+" : "-", ", ", norm(G, i), ")")
      if i < length(G)
        print(io, " ")
      end
    end
  else
    for i in 1:length(G)
      print(io, "(", scale(G, i), ", ", rank(G, i),
            ", ", det(G, i) == 1 ? "+" : "-",  ")")
      if i < length(G)
        print(io, " ")
      end
    end
  end
end

################################################################################
#
#  Basic properties
#
################################################################################

@doc Markdown.doc"""
    scale(G::LocalGenusHerm, i::Int) -> Int

Given a genus symbol for Hermitian lattices over $E/K$, return the $\mathfrak
P$-valuation of the scale of the $i$th Jordan block of $G$, where $\mathfrak
P$ is a prime ideal lying over $\mathfrak p$.
"""
scale(G::LocalGenusHerm, i::Int) = G.data[i][1]

@doc Markdown.doc"""
    scales(G::LocalGenusHerm, i::Int) -> Vector{Int}

Given a genus symbol for Hermitian lattices over $E/K$, return the $\mathfrak
P$-valuation of the scales of the Jordan blocks of $G$, where $\mathfrak
P$ is a prime ideal lying over $\mathfrak p$.
"""
scales(G::LocalGenusHerm) = map(i -> scale(G, i), 1:length(G))::Vector{Int}

@doc Markdown.doc"""
    rank(G::LocalGenusHerm, i::Int)

Given a genus symbol for Hermitian lattices over $E/K$, return the rank of the
$i$th Jordan block of $G$.
"""
rank(G::LocalGenusHerm, i::Int) = G.data[i][2]

@doc Markdown.doc"""
    rank(G::LocalGenusHerm) -> Int

Given a genus symbol $G$, return the rank of a lattice in $G$.
"""
function rank(G::LocalGenusHerm)
  return reduce(+, (rank(G, i) for i in 1:length(G)), init = Int(0))
end

@doc Markdown.doc"""
    ranks(G::LocalGenusHerm)

Given a genus symbol for Hermitian lattices over $E/K$, return the ranks of the
Jordan blocks of $G$.
"""
ranks(G::LocalGenusHerm) = map(i -> rank(G, i), 1:length(G))::Vector{Int}

@doc Markdown.doc"""
    det(G::LocalGenusHerm, i::Int) -> Int

Given a genus symbol for Hermitian lattices over $E/K$, return the determinant of the
$i$th Jordan block of $G$. This will be `1` or `-1` depending on whether the
determinant is a local norm or not.
"""
det(G::LocalGenusHerm, i::Int) = G.data[i][3]

@doc Markdown.doc"""
    discriminant(G::LocalGenusHerm, i::Int) -> Int

Given a genus symbol for Hermitian lattices over $E/K$, return the discriminant
of the $i$th Jordan block of $G$. This will be `1` or `-1` depending on whether
the discriminant is a local norm or not.
"""
function discriminant(G::LocalGenusHerm, i::Int)
  d = det(G)
  r = rank(G, i) % 4
  if r == 0 || r == 1
    return d
  end
  E = base_field(G)
  fl = islocal_norm(E, base_field(E)(-1), prime(G))
  if fl
    return d
  else
    return -d
  end
end

@doc Markdown.doc"""
    discriminant(G::LocalGenusHerm) -> Int

Given a genus symbol $G$, return the discriminant of a lattice in $G$. This will be
`1` or `-1` depending on whether the discriminant is a local norm or not.
"""
function discriminant(G::LocalGenusHerm)
  d = det(G)
  r = rank(G) % 4
  if r == 0 || r == 1
    return d
  end
  E = base_field(G)
  fl = islocal_norm(E, base_field(K)(-1), prime(G))
  if fl
    return d
  else
    return -d
  end
end

@doc Markdown.doc"""
    det(G::LocalGenusHerm) -> Int

Given a genus symbol $G$, return the determinant of a lattice in $G$. This will be
`1` or `-1` depending on whether the determinant is a local norm or not.
"""
function det(G::LocalGenusHerm)
  return reduce(*, (det(G, i) for i in 1:length(G)), init = Int(1))
end

@doc Markdown.doc"""
    dets(G::LocalGenusHerm) -> Vector{Int}

Given a genus symbol for Hermitian lattices over $E/K$, return the determinants
of the Jordan blocks of $G$. These will be `1` or `-1` depending on whether the
determinant is a local norm or not.
"""
dets(G::LocalGenusHerm) = map(i -> det(G, i), 1:length(G))::Vector{Int}

@doc Markdown.doc"""
    norm(G::LocalGenusHerm, i::Int) -> Int

Given a dyadic genus symbol for Hermitian lattices over $E/K$ at a prime
$\mathfrak p$, return the $\mathfrak p$-valuation of the norm of the $i$th Jordan
block of $G$.
"""
norm(G::LocalGenusHerm, i::Int) = begin @assert isdyadic(G); G.norm_val[i] end # this only works if it is dyadic

#@doc Markdown.doc"""
#    norms(G::LocalGenusHerm) -> Vector{Int}
#
#Given a dyadic genus symbol for Hermitian lattices over $E/K$ at a prime
#$\mathfrak p$, return the $\mathfrak p$-valuations of the norms of the Jordan
#blocks of $G$.
#"""

@doc Markdown.doc"""
    isramified(G::LocalGenusHerm) -> Bool

Given a genus symbol for Hermitian lattices at a prime $\mathfrak p$, return
whether $\mathfrak p$ is ramified.
"""
isramified(G::LocalGenusHerm) = G.isramified

@doc Markdown.doc"""
    isdyadic(G::LocalGenusHerm) -> Bool

Given a genus symbol for Hermitian lattices at a prime $\mathfrak p$, return
whether $\mathfrak p$ is dyadic.
"""
isdyadic(G::LocalGenusHerm) = G.isdyadic

#data(G::LocalGenusHerm) = G.data

@doc Markdown.doc"""
    length(G::LocalGenusHerm) -> Int

Return the number of Jordan blocks in a Jordan decomposition of $G$.
"""
length(G::LocalGenusHerm) = length(G.data)

@doc Markdown.doc"""
    base_field(G::LocalGenusHerm) -> NumField

Given a local genus symbol of Hermitian lattices over $E/K$, return $E$.
"""
base_field(G::LocalGenusHerm) = G.E

@doc Markdown.doc"""
    prime(G::LocalGenusHerm) -> NfOrdIdl

Given a local genus symbol of Hermitian lattices at a prime $\mathfrak p$,
return $\mathfrak p$.
"""
prime(G::LocalGenusHerm) = G.p

################################################################################
#
#  A local non-norm
#
################################################################################

function _non_norm_rep(G::LocalGenusHerm)
  if isdefined(G, :non_norm_rep)
    return G.non_norm_rep::elem_type(base_field(base_field(G)))
  end

  z = _non_norm_rep(base_field(G), base_field(base_field(G)), prime(G))
  G.non_norm_rep = z
  return z::elem_type(base_field(base_field(G)))
end

################################################################################
#
#  A local unifomizer
#
################################################################################

@doc Markdown.doc"""
    uniformizer(G::LocalGenusHerm) -> NumFieldElem

Given a local genus symbol of Hermitian lattices over $E/K$ at a prime $\mathfrak p$,
return a generator for the largest invariant ideal of $E$ containing $\mathfrak p$.
"""
function uniformizer(G::LocalGenusHerm)
  E = base_field(G)
  K = base_field(E)
  if isramified(G)
    lP = prime_decomposition(maximal_order(E), prime(G))
    @assert length(lP) == 1 && lP[1][2] == 2
    Q = lP[1][1]
    pi = p_uniformizer(Q)
    A = automorphisms(E)
    uni = A[1](elem_in_nf(pi)) * A[2](elem_in_nf(pi))
    @assert iszero(coeff(uni, 1))
    @assert islocal_norm(E, coeff(uni , 0), prime(G))
    return coeff(uni, 0)
  else
    return K(uniformizer(prime(G)))
  end
end

################################################################################
#
#  Additional dyadic information
#
################################################################################

# Get the "ni" for the ramified dyadic case
function _get_ni_from_genus(G::LocalGenusHerm)
  @assert isramified(G)
  @assert isdyadic(G)
  t = length(G)
  z = Vector{Int}(undef, t)
  for i in 1:t
    ni = minimum(2 * max(0, scale(G, i) - scale(G, j)) + 2 * norm(G, j) for j in 1:t)
    z[i] = ni
  end
  return z
end

################################################################################
#
#  Determinant representative
#
################################################################################

# TODO: I don't understand what this is doing. I know that d tells me if it is
#       a norm or not. Why do I have to multiply with a uniformizer
@doc Markdown.doc"""
    det_representative(G::LocalGenusHerm) -> NumFieldElem

Return a representative for the norm class of the determinant of $G$.
"""
function det_representative(G::LocalGenusHerm)
  z = G.data
  d = det(G)
  v = sum(scale(G, i) * rank(G, i) for i in 1:length(G))
  if isramified(G)
    v = div(v, 2)
  end
  if d == 1
    u = one(base_field(base_field(G)))
  else
    @assert isramified(G)
    u = _non_norm_rep(G)
  end
  return u * uniformizer(G)^v
end

################################################################################
#
#  Gram matrix
#
################################################################################

@doc Markdown.doc"""
    gram_matrix(G::LocalGenusHerm)

Return a matrix $M$, such that a lattice with Gram matrix $M$ is an element of
the given genus.
"""
function gram_matrix(G::LocalGenusHerm)
  return diagonal_matrix(dense_matrix_type(base_field(G))[gram_matrix(G, i) for i in 1:length(G)])
end

function gram_matrix(G::LocalGenusHerm, l::Int)
  i = scale(G, l)
  m = rank(G, l)
  d = det(G, l)
  E = base_field(G)
  K = base_field(E)
  p = elem_in_nf(p_uniformizer(prime(G)))
  A = automorphisms(E)
  _a = gen(E)
  conj = A[1](_a) == _a ? A[2] : A[1]

  if !isramified(G)
    return diagonal_matrix([E(p)^i for j in 1:m])
  end

  # ramified

  lQ = prime_decomposition(maximal_order(E), prime(G))
  @assert length(lQ) == 1 && lQ[1][2] == 2
  Q = lQ[1][1]
  pi = elem_in_nf(p_uniformizer(Q))
  H = matrix(E, 2, 2, elem_type(E)[zero(E), pi^i, conj(pi)^i, zero(E)])

  if !isdyadic(G)
    # non-dyadic
    if iseven(i)
      # According to Kir16, there the last exponent should be i/2 * (1 - m)
      if d == 1
        u = one(K)
      else
        u = _non_norm_rep(G)
      end
      if m == 1
        fl = islocal_norm(E, p^div(i, 2), prime(G))
        if d == -1
          return diagonal_matrix([(fl ? _non_norm_rep(G) : one(K)) * E(p)^div(i, 2)])
        else
          return diagonal_matrix([(!fl ? _non_norm_rep(G) : one(K)) * E(p)^div(i, 2)])
        end
      end
      return diagonal_matrix(push!(elem_type(E)[E(p)^div(i, 2) for j in 1:(m - 1)], u * E(p)^(div(i, 2) * (m - 1))))
    else
      return diagonal_matrix(dense_matrix_type(E)[H for j in 1:div(m, 2)])
    end
  end

  # dyadic

  k = norm(G, l)

  # I should cache this e
  e = valuation(different(maximal_order(E)), Q)

  if isodd(m)
    # odd rank
    @assert 2*k == i
    r = div(m - 1, 2)
    nn = mod(m, 4)
    if d == 1
      discG = K(1)
    else
      discG = elem_in_nf(_non_norm_rep(G))
    end
    if nn == 0 || nn == 1
      discG = discG
    else
      discG = -discG
    end
    @assert iseven(i)
    if islocal_norm(E, discG*p^(-m * div(i, 2)), prime(G))
      u = K(1)
    else
      u = _non_norm_rep(G)
    end

    U = matrix(E, 1, 1, [u * p^(div(i, 2))])
    return diagonal_matrix(append!(typeof(U)[U], [H for j in 1:r]))
  else
    # even rank
    r = div(m, 2) - 1
    if islocal_norm(E, K((-1)^div(m, 2)), prime(G)) == (d == 1)
      # hyperbolic
      @assert i + e >= 2 * k
      @assert 2 * k >= i
      U = matrix(E, 2, 2, [p^k, pi^i, conj(pi)^i, 0])
      return diagonal_matrix(append!(typeof(U)[U], [H for j in 1:r]))
    else # not hyperbolic
      @assert i + e > 2 * k
      @assert 2 * k >= i
      u = _non_norm_rep(G)
      u0 = u - 1
      U = matrix(E, 2, 2, [p^k, pi^i, conj(pi)^i, -p^(i -k) * u0])
      return diagonal_matrix(append!(typeof(U)[U], [H for j in 1:r]))
    end
  end
end

################################################################################
#
#  Representative
#
################################################################################

@doc Markdown.doc"""
    representative(G::LocalGenusHerm) -> HermLat

Given a local genus, return a Hermitian lattice contained in this genus.
"""
function representative(G::LocalGenusHerm)
  E = G.E
  return lattice(hermitian_space(E, gram_matrix(G)))
end

################################################################################
#
#  Constructor
#
################################################################################

@doc Markdown.doc"""
    genus(HermLat, E::NumField, p::Idl, data::Vector{Tuple{Int, Int, Int}};
                                        type = :det)
                                                              -> LocalGenusHerm

Construct the local genus symbol of hermitian lattices over $E$ at the prime ideal
$\mathfrak p$ with the invariants specified by `data`.

If the prime ideal is good, the vector `data` contains for each block of the
Jordan decomposition a pair `(s, r, d)`, where `s` is the scale, `r` the
rank. The value `d` must be in `[-1, 1]` and indicates whether the determinant
of the block is a local norm or not.

If the optional `type` keyword is set to `:disc`, then `d` is interpreted as the
norm class of the discriminant of the corresponding Jordan block.
"""
genus(::Type{HermLat}, E, p, data; type)

# rank zero genus
function genus(::Type{HermLat}, E::S, p::T) where {S, T}
  if isramified(maximal_order(E), p) && isdyadic(p)
    return genus(HermLat, E, p, Vector{Tuple{Int, Int, Int, Int}}())
  else
    return genus(HermLat, E, p, Vector{Tuple{Int, Int, Int}}())
  end
end

function genus(::Type{HermLat}, E::S, p::T, data::Vector{Tuple{Int, Int, Int}}; type = :det) where {S <: NumField, T}
  z = LocalGenusHerm{S, T}()
  z.E = E
  z.p = p
  z.isdyadic = isdyadic(p)
  z.isramified = isramified(maximal_order(E), p)
  @assert !(isramified(z) && isdyadic(z))

  if type !== :det && type !== :disc
    throw(error("type :$type must be :disc or :det"))
  end

  if !z.isramified || type === :det
    z.data = copy(data)
  else
    type !== :disc && throw(error("type :$type must be :disc or :det"))
    fl = islocal_norm(E, base_field(E)(-1), p)
    if fl
      z.data = copy(data)
    end
    # Now -1 is not a local norm, so we adjust whenever the rank is 2, 3 mod 4.
    z.data = Vector{Tuple{Int, Int, Int}}(undef, length(data))
    for i in 1:length(data)
      r = data[i][2] % 4
      if r == 0 || r == 1
        z.data[i] = data[i]
      else
        z.data[i] = (data[i][1], data[i][2], (-1) * data[i][3])
      end
    end
  end
  return z
end

function genus(::Type{HermLat}, E::S, p::T, data::Vector{Tuple{Int, Int, Int, Int}}; type = :det) where {S <: NumField, T}
  z = LocalGenusHerm{S, T}()
  z.E = E
  z.p = p
  z.isdyadic = isdyadic(p)
  z.isramified = isramified(maximal_order(E), p)
  if isramified(z) && isdyadic(z)
    z.data = Tuple{Int, Int, Int}[Base.front(v) for v in data]
    z.norm_val = Int[v[end] for v in data]
    z.ni = _get_ni_from_genus(z)

    for i in 1:length(z.data)
      # If the rank is odd, then n(L) * O_E = s(L), so n = 2 * s,
      # since n is the valuation in K and the extension is ramified.
      v = z.data[i]
      if isodd(v[2])
        if 2 * z.norm_val[i] != v[1]
          throw(error("""Not a valid local genus in block $(i):
                      Scale ($(v[1])) must be twice the norm ($(z.norm_val[i]))"""))
        end
      end
    end
  else
    z.data = Tuple{Int, Int, Int}[Base.front(v) for v in data]
  end

  if type !== :det && type !== :disc
    throw(error("type :$type must be :disc or :det"))
  end

  if type === :disc
    fl = islocal_norm(E, base_field(E)(-1), p)
    if !fl
      for i in 1:length(z.data)
        r = z.data[i][2] % 4
        if r == 0 || r == 1
          continue
        else
          z.data[i] = (z.data[i][1], z.data[i][2], (-1) * z.data[i][3])
        end
      end
    end
  end

  return z::LocalGenusHerm{S, T}
end

################################################################################
#
#  Genus symbol of a lattice
#
################################################################################

# TODO: better documentation

@doc Markdown.doc"""
    genus(L::HermLat, p::NfOrdIdl) -> LocalGenusHerm

Returns the genus of $L$ at the prime ideal $\mathfrak p$.

See [Kir16, Definition 8.3.1].
"""
function genus(L::HermLat, q)
  c = get_special(L, :local_genera)
  S = ideal_type(base_ring(base_ring(L)))
  if c === nothing
    symbols = Dict{S, LocalGenusHerm{typeof(base_field(L)), S}}()
    set_special(L, :local_genera => symbols)
  else
    symbols = c::Dict{S, LocalGenusHerm{typeof(base_field(L)), S}}
    if order(q) !== base_ring(base_ring(L))
      if haskey(symbols, minimum(q))
        return symbols[minimum(q)]
      end
    elseif q isa fmpz
      qq = ideal(base_ring(base_ring(L)), q)
      if haskey(symbols, qq)
        return symbols[qq]
      end
    else
      if haskey(symbols, q)
        return symbols[q]
      end
    end
  end

  if order(q) !== base_ring(base_ring(L))
    g = _genus(L, minimum(q))
    symbols[minimum(q)] = g
  elseif q isa fmpz
    qq = ideal(base_ring(base_ring(L)), q)
    g = _genus(L, qq)
    symbols[qq] = g
    return g
  else
    g = _genus(L, q)
    symbols[q] = g
    return g
  end
end

function _genus(L::HermLat, p)
  sym = _genus_symbol(L, p)
  G = genus(HermLat, nf(base_ring(L)), p, sym)
  # Just for debugging
  @hassert :Lattice 1 begin
    if isdyadic(G) && isramified(G)
      GG = _genus_symbol_kirschmer(L, p)
      all(i -> GG[i][4] == G.ni[i], 1:length(G))
    else
      true
    end
  end
  return G
end

################################################################################
#
#  Test if lattice is contained in local genus
#
################################################################################

@doc Markdown.doc"""
    in(L::HermLat, G::LocalGenusHerm) -> Bool

Test if the lattice $L$ is contained in the local genus $G$.
"""
Base.in(L::HermLat, G::LocalGenusHerm) = genus(L, prime(G)) == G

################################################################################
#
#  Equality
#
################################################################################

@doc Markdown.doc"""
    ==(G1::LocalGenusHerm, G2::LocalGenusHerm)

Test if two local genera are equal.
"""
function ==(G1::LocalGenusHerm, G2::LocalGenusHerm)
  if base_field(G1) != base_field(G2)
    return false
  end

  if prime(G1) != prime(G2)
    return false
  end

  if length(G1) != length(G2)
    return false
  end

  t = length(G1)

  p = prime(G1)

  # We now check the Jordan type

  if any(i -> scale(G1, i) != scale(G2, i), 1:t)
    return false
  end

  if any(i -> (rank(G1, i) != rank(G2, i)), 1:t)
    return false
  end

  if det(G1) != det(G2) # rational spans must be isometric
    return false
  end

  if any(i -> (rank(G1, i) != rank(G2, i)), 1:t)
    return false
  end

  if !isramified(G1) # split or unramified
    return true
    # Everything is normal and the Jordan decomposition types agree
    #return all(det(G1, i) == det(G2, i) for i in 1:t)
  end

  if isramified(G1) && !isdyadic(G1) # ramified, but not dyadic
    # If s(L_i) is odd, then L_i = H(s(L_i))^(rk(L_i)/2) = H(s(L_i'))^(rk(L_i')/2) = L_i'
    # So all L_i, L_i' are locally isometric, in particular L_i is normal if and only if L_i' is normal
    # If s(L_i) = s(L_i') is even, then both L_i and L_i' have orthgonal bases, so they are
    # in particular normal.

    # Thus we only have to check Theorem 3.3.6 4.
    return all(i -> det(G1, i) == det(G2, i), 1:t)
    # TODO: If think I only have to do something if the scale is even. Check this!
  end

  # Dyadic ramified case

  # First check if L_i is normal if and only if L_i' is normal, i.e.,
  # that the Jordan decompositions agree
  for i in 1:t
    if (scale(G1, i) == 2 * norm(G1, i)) != (scale(G2, i) == 2 * norm(G2, i))
      return false
    end
  end

  if any(i -> G1.ni[i] != G2.ni[i], 1:t)
    return false
  end

  E = base_field(G1)
  lQ = prime_decomposition(maximal_order(E), p)
  @assert length(lQ) == 1 && lQ[1][2] == 2
  Q = lQ[1][1]

  e = valuation(different(maximal_order(E)), Q)

  for i in 1:(t - 1)
    dL1prod = prod(det(G1, j) for j in 1:i)
    dL2prod = prod(det(G2, j) for j in 1:i)

    d = dL1prod * dL2prod

    if d != 1
      if 2 * (e - 1) < G1.ni[i] + G1.ni[i + 1] - 2 * scale(G1, i)
        return false
      end
    end
  end

  return true
end

function _genus_symbol(L::HermLat, q)
  if order(q) != base_ring(base_ring(L))
    return __genus_symbol(L, minimum(q))
  else
    return __genus_symbol(L, q)
  end
  end

function __genus_symbol(L::HermLat, p)
  @assert order(p) == base_ring(base_ring(L))
  B, G, S = jordan_decomposition(L, p)
  R = base_ring(L)
  E = nf(R)
  K = base_field(E)
  local sym::Vector{Tuple{Int, Int, Int, Int}}
  if !isdyadic(p) || !isramified(R, p)
    # The last entry is a dummy to make the compiler happier
    sym = Tuple{Int, Int, Int, Int}[ (S[i], nrows(B[i]), islocal_norm(E, coeff(det(G[i]), 0), p) ? 1 : -1, 0) for i in 1:length(B)]
  else
    P = prime_decomposition(R, p)[1][1]
    pi = E(K(uniformizer(p)))
    sym = Tuple{Int, Int, Int, Int}[]
    for i in 1:length(B)
      normal = _get_norm_valuation_from_gram_matrix(G[i], P) == S[i]
     GG = diagonal_matrix(dense_matrix_type(E)[pi^(max(0, S[i] - S[j])) * G[j] for j in 1:length(B)])
      #v = _get_norm_valuation_from_gram_matrix(GG, P)
      #@assert v == valuation(R * norm(lattice(hermitian_space(E, GG), identity_matrix(E, nrows(GG)))), P)
      r = nrows(B[i]) # rank
      s = S[i] # P-valuation of scale(L_i)
      det_class = islocal_norm(E, coeff(det(G[i]), 0), p) ? 1 : -1  # Norm class of determinant
      normi = _get_norm_valuation_from_gram_matrix(G[i], P) # P-valuation of norm(L_i)
      @assert mod(normi, 2) == 0 # I only want p-valuation
      push!(sym, (s, r, det_class, div(normi, 2)))
    end
  end
  return sym
end

################################################################################
#
#  Global genus
#
################################################################################

mutable struct GenusHerm{S, T, U, V}
  E::S
  primes::Vector{T}
  LGS::Vector{U}
  rank::Int
  signatures::V

  function GenusHerm(E::S, r, LGS::Vector{U}, signatures::V) where {S, U, V}
    K = base_field(E)
    primes = Vector{ideal_type(order_type(K))}(undef, length(LGS))

    for i in 1:length(LGS)
      primes[i] = prime(LGS[i])
      @assert r == rank(LGS[i])
    end
    z = new{S, eltype(primes), U, V}(E, primes, LGS, r, signatures)
    return z
  end
end

genus_herm_type(E) = GenusHerm{typeof(E), ideal_type(order_type(base_field(E))), local_genus_herm_type(E), Dict{InfPlc, Int}}

@doc Markdown.doc"""
    ==(G1::GenusHerm, G2::GenusHerm)

Test if two genera are equal.
"""
function ==(G1::GenusHerm, G2::GenusHerm)
  if G1.E != G2.E
    return false
  end

  if length(G1.primes) != length(G2.primes)
    return false
  end

  if G1.signatures != G2.signatures
    return false
  end

  for p in G1.primes
    if !(p in G2.primes)
      return false
    end
    if G1[p] != G2[p]
      return false
    end
  end

  return true
end

function Base.show(io::IO, G::GenusHerm)
  print(io, "Genus symbol over")
  print(io, G.E)
  print(io, "\n", "and local genera",)
  for g in G.LGS
    print(io, "\n")
    print(IOContext(io, :compact => true), prime(g), " => ", g)
  end
  print(io, "\n and signature")
  for (pl, v) in G.signatures
    print(io, "\n")
    _print_short(io, pl.r)
    print(io, " => ")
    print(io, v)
  end
end

function _print_short(io::IO, a::arb)
  r = BigFloat(a)
  s = string(r)
  if length(s) >= 10
    ss = s[1:9] * "…"
  else
    ss = s
  end
  print(io, ss)
end

function _print_short(io::IO, a::acb)
  _print_short(io, real(a))
  if !iszero(imag(a))
    print(io, " + ")
    _print_short(io, imag(a))
    print(io, " * i")
  end
end

################################################################################
#
#  Sum of genera
#
################################################################################

function orthogonal_sum(G1::GenusHerm, G2::GenusHerm)
  @req G1.E === G2.E "Genera must have same base field"
  E = G1.E
  LGS = local_genus_herm_type(G1.E)[]
  prim = eltype(primes(G1))[]
  P1 = Set(primes(G1))
  P2 = Set(primes(G2))
  for p in union(P1, P2)
    if p in P1
      i = findfirst(g -> prime(g) == p, G1.LGS)
      g1 = G1.LGS[i]
    else
      @assert !isramified(maximal_order(E), p)
      g1 = genus(HermLat, p, [(0, rank(G1), 1)])
    end

    if p in P2
      i = findfirst(g -> prime(g) == p, G2.LGS)
      g2 = G2.LGS[i]
    else
      @assert !isramified(maximal_order(E), p)
      g2 = genus(HermLat, p, [(0, rank(G2), 1)])
    end

    g3 = orthogonal_sum(g1, g2)
    push!(prim, p)
    push!(LGS, g3)
  end
  sig1 = G1.signatures
  sig2 = G2.signatures
  g3 = merge(+, sig1, sig2)
  return genus(LGS, g3)
end

Base.:(+)(G1::GenusHerm, G2::GenusHerm) = orthogonal_sum(G1, G2)

################################################################################
#
#  Test if lattice is contained in genus
#
################################################################################

@doc Markdown.doc"""
    in(L::HermLat, G::GenusHerm) -> Bool

Test if the lattice $L$ is contained in the genus $G$.
"""
Base.in(L::HermLat, G::GenusHerm) = genus(L) == G

# This could be made more efficient

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

function Base.show(io::IO, ::MIME"text/plain", G::GenusHerm)
  print(io, "Global genus symbol\n")
  for i in 1:length(G.primes)
    print(IOContext(io, :compact => true), G.primes[i], " => ", G.LGS[i],)
    if i < length(G.primes)
      print(io, "\n")
    end
  end
  for (pl, v) in G.signatures
    print(io, "\n")
    _print_short(io, pl.r)
    print(io, " => ")
    print(io, v)
  end
end

################################################################################
#
#  Constructor
#
################################################################################

@doc Markdown.doc"""
    genus(L::Vector{LocalGenusHerm},
          signatures::Vector{Tuple{InfPlc, Int}}) -> GenusHerm

Given a set of local genera and signatures, return the corresponding global
genus.
"""
function genus(L::Vector{<:LocalGenusHerm}, signatures::Dict{InfPlc, Int})
  @assert !isempty(L)
  @assert all(N >= 0 for (_, N) in signatures)
  if !_check_global_genus(L, signatures)
    throw(error("Invariants violate the product formula."))
  end
  r = rank(first(L))
  E = base_field(first(L))
  return GenusHerm(E, r, L, signatures)
end

function genus(L::Vector, signatures::Vector{Tuple{InfPlc, Int}})
  return genus(L, Dict(signatures))
end

@doc Markdown.doc"""
    genus(L::HermLat) -> GenusHerm

Given a Hermitian lattice, return the genus it belongs to.
"""
function genus(L::HermLat)
  c = get_special(L, :genus)
  S = ideal_type(base_ring(base_ring(L)))
  if c === nothing
    G = _genus(L)
    set_special(L, :genus => G)
    return G
  else
    return c::genus_herm_type(base_field(L))
  end
end

function _genus(L::HermLat)
  bad = bad_primes(L)
  for p in support(discriminant(base_ring(L)))
    if !(p in bad)
      push!(bad, p)
    end
  end

  for p in support(2 * base_ring(base_ring(L)))
    if !(p in bad)
      push!(bad, p)
    end
  end

  S = real_places(base_field(base_field(L)))
  D = diagonal(rational_span(L))
  signatures = Dict{InfPlc, Int}(s => count(d -> isnegative(d, s), D) for s in S)
  return genus([genus(L, p) for p in bad], signatures)
end

################################################################################
#
#  Basic access
#
################################################################################

@doc Markdown.doc"""
    primes(G::GenusHerm) -> Vector{NfOrdIdl}

Return the primes of a global genus symbol.
"""
primes(G::GenusHerm) = G.primes

@doc Markdown.doc"""
    signatures(G::GenusHerm) -> Dict{InfPlc, Int}

Return the signatures at the infinite places, which for each real place is
given by the negative index of inertia.
"""
signatures(G::GenusHerm) = G.signatures

################################################################################
#
#  Consistency check
#
################################################################################

function _check_global_genus(LGS, signatures)
  _non_norm = _non_norm_primes(LGS)
  P = length(_non_norm)
  I = length([(s, N) for (s, N) in signatures if mod(N, 2) == 1])
  if mod(P + I, 2) == 1
    return false
  end
  return true
end

################################################################################
#
#  Primes at which the determinant is not a local norm
#
################################################################################

function _non_norm_primes(LGS::Vector{LocalGenusHerm{S, T}}) where {S, T}
  z = T[]
  for g in LGS
    p = prime(g)
    d = det(g)
    if d != 1
      push!(z, p)
    end
  end
  return z
end

################################################################################
#
#  Convenient functions
#
################################################################################

function Base.getindex(G::GenusHerm, P)
  i = findfirst(isequal(P), G.primes)
  if i === nothing
    throw(error("No local genus symbol at $P"))
  end
  return G.LGS[i]
end

################################################################################
#
#  Compute representatives of genera
#
################################################################################

function _hermitian_form_with_invariants(E, dim, P, N)
  #@show dim, P, N
  #@show E, dim, N, [minimum(p) for p in P]
  K = base_field(E)
  R = maximal_order(K)
#  require forall{n: n in N | n in {0..dim}}: "Number of negative entries is impossible";
  infinite_pl = [ p for p in real_places(K) if length(infinite_places(E, p)) == 1 ]
  length(N) != length(infinite_pl) && error("Wrong number of real places")
  S = maximal_order(E)
  prim = [ p for p in P if length(prime_decomposition(S, p)) == 1 ] # only take non-split primes
  #@show prim
  I = [ p for p in keys(N) if isodd(N[p]) ]
  !iseven(length(I) + length(P)) && error("Invariants do not satisfy the product formula")
  e = gen(E)
  x = 2 * e - trace(e)
  b = coeff(x^2, 0) # b = K(x^2)
  #@show b, prim, I
  a = _find_quaternion_algebra(b, prim, I)
  #@show a
  D = elem_type(E)[]
  for i in 1:(dim - 1)
    if length(I) == 0
      push!(D, one(E))
    else
      el = E(_weak_approximation(infinite_pl, [N[p] >= i ? -1 : 1 for p in infinite_pl]))
      push!(D, el)
    end
  end
  push!(D, a * prod(D))
  Dmat = diagonal_matrix(D)
  #@show Dmat
  dim0, P0, N0 = _hermitian_form_invariants(Dmat)
  #@show P0
  #@show prim
  @assert dim == dim0
  @assert Set(prim) == Set(P0)
  @assert N == N0
  return Dmat
end

function _hermitian_form_invariants(M)
  E = base_ring(M)
  K = base_field(E)
  @assert degree(E) == 2
  A = automorphisms(E)
  a = gen(E)
  v = A[1](a) == a ? A[2] : A[1]

  @assert M == transpose(_map(M, v))
  d = coeff(det(M), 0) # K(det(M))
  P = Dict()
  for p in keys(factor(d * maximal_order(K)))
    if islocal_norm(E, d, p)
      continue
    end
    P[p] = true
  end
  for p in keys(factor(discriminant(maximal_order(E))))
    if islocal_norm(E, d, p)
      continue
    end
    P[p] = true
  end
  D = diagonal(_gram_schmidt(M, v)[1])
  I = Dict([ p=>length([coeff(d, 0) for d in D if isnegative(coeff(d, 0), p)]) for p in real_places(K) if length(infinite_places(E, p)) == 1])
  return ncols(M), collect(keys(P)), I
end

base_field(G::GenusHerm) = G.E

rank(G::GenusHerm) = G.rank

function representative(G::GenusHerm)
  P = _non_norm_primes(G.LGS)
  E = base_field(G)
  V = hermitian_space(E, _hermitian_form_with_invariants(base_field(G), rank(G), P, G.signatures))
  @vprint :Lattice 1 "Finding maximal integral lattice\n"
  M = maximal_integral_lattice(V)
  lp = G.primes
  for g in G.LGS
    p = prime(g)
    @vprint :Lattice 1 "Finding representative for $g at $(prime(g))...\n"
    L = representative(g)
    @hassert :Lattice 1 genus(L, p) == g
    #@show coefficient_ideals(pseudo_matrix(L))
    #@show matrix(pseudo_matrix(L))
    @vprint :Lattice 1 "Finding sublattice\n"
    M = locally_isometric_sublattice(M, L, p)
  end
  return M
end

################################################################################
#
#  Enumeration of local genera
#
################################################################################

@doc Markdown.doc"""
    local_genera_hermitian(E::NumField, p::NfOrdIdl, rank::Int,
                 det_val::Int, max_scale::Int) -> Vector{LocalGenusHerm}

Return all local genera of Hermitian lattices over $E$ at $\mathfrak p$ with
rank `rank`, scale valuation bounded by `max_scale` and determinant valuation
equal to `det_val`.
"""
function local_genera_hermitian(E, p, rank::Int, det_val::Int, max_scale::Int, is_ramified = isramified(maximal_order(E), p))
  #@show E, p, rank, det_val, max_scale, is_ramified
  is_inert = !is_ramified && length(prime_decomposition(maximal_order(E), p)) == 1
  if is_ramified
    # the valuation is with respect to p
    # but the scale is with respect to P
    # in the ramified case p = P**2 and thus
    # the determinant of a block is
    # P^(scale*rank) = p^(scale*rank/2)
    det_val *= 2
  end

  K = number_field(order(p))

  scales_rks = Vector{Tuple{Int, Int}}[] # possible scales and ranks

  for rkseq in _integer_lists(rank, max_scale + 1)
    d = 0
    pgensymbol = Tuple{Int, Int}[]
    for i in 0:(max_scale + 1) - 1
      d += i * rkseq[i + 1]
      if rkseq[i + 1] != 0
        push!(pgensymbol, (i, rkseq[i + 1]))
      end
    end
    if d == det_val
        push!(scales_rks, pgensymbol)
    end
  end

  if !is_ramified
    # I add the 0 to make the compiler happy
    symbols = Vector{LocalGenusHerm{typeof(E), typeof(p)}}(undef, length(scales_rks))
    for i in 1:length(scales_rks)
      g = scales_rks[i]
      z = Tuple{Int, Int, Int, Int}[]
      for b in g
        # We have to be careful.
        # If p is inert, then the norm is not surjective.
        if !is_inert || iseven(b[1] * b[2])
          push!(z, (b[1], b[2], 1, 0))
        else
          push!(z, (b[1], b[2], -1, 0))
        end
      end
      symbols[i] = genus(HermLat, E, p, z)
    end
    return symbols
  end

  scales_rks = Vector{Tuple{Int, Int}}[g for g in scales_rks if all((mod(b[1]*b[2], 2) == 0) for b in g)]

  symbols = LocalGenusHerm{typeof(E), typeof(p)}[]
  #hyperbolic_det = hilbert_symbol(K(-1), gen(K)^2//4 - 1, p)
  hyperbolic_det = islocal_norm(E, K(-1), p) ? 1 : -1
  if !isdyadic(p) # non-dyadic
    for g in scales_rks
      n = length(g)
      dets = Vector{Int}[]
      for b in g
        if mod(b[1], 2) == 0
          push!(dets, [1, -1])
        end
        if mod(b[1], 2) == 1
          push!(dets, [hyperbolic_det^(div(b[2], 2))])
        end
      end

      for d in cartesian_product_iterator(dets, inplace = true)# Iterators.product(dets...)
        g2 = Vector{Tuple{Int, Int, Int, Int}}(undef, length(g))
        for k in 1:n
          # Again the 0 for dummy purposes
          g2[k] = (g[k]..., d[k], 0)
        end
        push!(symbols, genus(HermLat, E, p, g2))
      end
    end
    return symbols
  end

  # Ramified case
  lp = prime_decomposition(maximal_order(E), p)
  @assert length(lp) == 1 && lp[1][2] == 2
  P = lp[1][1]

  e = valuation(different(maximal_order(E)), P)
  # only for debugging
  scales_rks = reverse(scales_rks)

  for g in scales_rks
    n = length(g)
    det_norms = Vector{Vector{Int}}[]
    for b in g
      if mod(b[2], 2) == 1
        @assert iseven(b[1])
        push!(det_norms, [[1, div(b[1], 2)], [-1, div(b[1], 2)]])
      end
      if mod(b[2], 2) == 0
        dn = Vector{Int}[]
        i = b[1]
        # (i + e) // 2 => k >= i/2
        k = ceil(Int, i//2)
        while 2 * k < i + e
          push!(dn, Int[1, k])
          push!(dn, Int[-1, k])
          k += 1
        end

        if mod(i + e, 2) == 0
          push!(dn, Int[hyperbolic_det^(div(b[2], 2)), k])
        end
        push!(det_norms, dn)
        #for k in (ceil(Int, Int(i)//2)):(div(Int(i + e), 2) - 1)
        #  push!(dn, Int[1, k])
        #  push!(dn, Int[-1, k])
        #end
        #push!(dn, Int[hyperbolic_det^(div(b[2], 2)), div(i + e, 2)])
        #if mod(i + e, 2) == 1
        #  push!(dn, Int[-hyperbolic_det^(div(b[2], 2)), div(i + e, 2)])
        #end
        #push!(det_norms, dn)
      end
    end

    for dn in cartesian_product_iterator(det_norms, inplace = true)
      g2 = Vector{Tuple{Int, Int, Int, Int}}(undef, length(g))
      for k in 1:n
        g2[k] = (g[k][1], g[k][2], dn[k][1], dn[k][2])
      end
      h = genus(HermLat, E, p, g2)
      if !(h in symbols)
        push!(symbols, h)
      end
    end
  end
  return symbols
end

@doc Markdown.doc"""
    genera_hermitian(E::NumField, rank::Int, determinant::Int, max_scale = nothing) -> Vector{GenusHerm}

Return all genera of Hermitian lattices over $E$ at $\mathfrak p$ with
rank `rank`, scale valuation bounded by `max_scale` and determinant valuation
equal to `det_val`.
"""
function genera_hermitian(E, rank, signatures, determinant; max_scale = nothing)
  K = base_field(E)
  OE = maximal_order(E)
  if max_scale === nothing
    _max_scale = determinant
  else
    _max_scale = max_scale
  end

  primes = support(discriminant(OE))
  for p in support(norm(determinant))
    if !(p in primes)
      push!(primes, p)
    end
  end

  local_symbols = Vector{local_genus_herm_type(E)}[]

  ms = norm(_max_scale)
  ds = norm(determinant)
  for p in primes
    det_val = valuation(ds, p)
    mscale_val = valuation(ms, p)
    det_val = div(det_val, 2)
    if !isramified(OE, p)
      mscale_val = div(mscale_val, 2)
    end
    push!(local_symbols, local_genera_hermitian(E, p, rank, det_val, mscale_val))
  end

  res = genus_herm_type(E)[]
  it = cartesian_product_iterator(local_symbols, inplace = true)
  for gs in it
    c = copy(gs)
    b = _check_global_genus(c, signatures)
    if b
      push!(res, GenusHerm(E, rank, c, signatures))
    end
  end

  return res
end

################################################################################
#
#  Genus representatives
#
################################################################################

# Return b, p, bad
# b = isdefinite(L)
# p = prime ideal of base_ring(L) which can be used for the neighbor method
# bad = bad primes of L, where L,p is not modular or p is dyadic and dividing discriminant(S)
function smallest_neighbour_prime(L)
  S = base_ring(L)
  R = base_ring(S)
  lp = bad_primes(L)
  bad = ideal_type(R)[p for p in lp if !ismodular(L, p)[1] ]
  for (p,_) in factor(discriminant(S))
    if isdyadic(p) && !(p in bad)
      push!(bad, p)
    end
  end

  if !isdefinite(L)
    return false, 1*S, bad
  end

  # TODO: This does not find the prime ideal with smallest norm,
  # but with smallest minimum ...

  m = rank(L)
  p = 1
  P = ideal_type(R)[]
  while length(P) == 0
    p = next_prime(p)
    lp = ideal_type(R)[ p[1] for p in prime_decomposition(R, p)]
    P = setdiff(lp, bad)
    if m == 2
      P = filter(p -> isisotropic(L, p), P)
    end
  end
  Q = prime_decomposition(S, P[1])[1][1]

  if isempty(bad)
    I = 1 * S
  else
    I = prod(bad) * S
  end
  n = absolute_norm(Q)
  if n >= 1000
    PP = prime_ideals_up_to(S, 1000)
    for QQ in PP
      if !iscoprime(QQ, I)
        continue
      end

      if isisotropic(L, QQ)
        return true, QQ, bad
      end
    end
  end
  PP = prime_ideals_up_to(S, n)
  for QQ in PP
    if !iscoprime(QQ, I)
      continue
    end
    if isisotropic(L, QQ)
      return true, QQ, bad
    end
  end
  throw(error("Impossible"))
end

function genus_generators(L::HermLat)
  R = base_ring(L)
  E = nf(R)
  D = different(R)
  b, P0, bad = smallest_neighbour_prime(L)

  local bad_prod::ideal_type(base_ring(R))

  if isempty(bad)
    bad_prod = 1 * base_ring(R)
  else
    bad_prod = prod(bad)
  end

  # First the ideals coming from the C/C0 quotient
  Eabs, EabstoE = absolute_simple_field(ambient_space(L))
  Rabs = maximal_order(Eabs)
  C, h = class_group(Rabs)
  RR = base_ring(R)
  C0 = support(D)
  CC, hh = class_group(RR)
  for p in find_gens(pseudo_inv(h), PrimesSet(2, -1))[1]
    if !(p in C0)
      push!(C0, sum(R * R(EabstoE(elem_in_nf(b))) for b in basis(p)))
    end
  end
  Q0, q0 = quo(C, elem_type(C)[ h\ideal(Rabs, [Rabs(EabstoE\b) for b in absolute_basis(i)]) for i in C0])
  q00 = pseudo_inv(q0) * h
  PP = ideal_type(R)[]

  local F::GaloisField

  local W::Generic.QuotientModule{gfp_elem}

  if iseven(rank(L))
    for (P, e) in factor(D)
      G = genus(L, P)
      if any(i -> isodd(rank(G, i)), 1:length(G))
        continue
      elseif !isdyadic(P)
        if any(i -> iseven(scale(G, i)), 1:length(G))
          continue
        end
      else
        if any(i -> isodd(e  + scale(G, i)) || (e + scale(G, i) != G.ni[i]), 1:length(G))
          continue
        end
      end
      push!(PP, P)
    end

    if !isempty(PP)
      U, f = unit_group_fac_elem(Rabs)
      UU, ff = unit_group_fac_elem(RR)
      nnorm = hom(U, UU, GrpAbFinGenElem[ff\FacElem(nf(RR)(norm(f(U[i])))) for i in 1:ngens(U)])
      l = length(PP)
      VD = Int[ valuation(D, P) for P in PP ]
      K, k = kernel(nnorm)
      F = GF(2, cached = false)
      V = VectorSpace(F, length(PP))
      S = elem_type(V)[]
      for u in gens(K)
        z = elem_type(F)[]
        for i in 1:length(PP)
          zz = R(EabstoE(evaluate(f(k(u))))) - 1
          if iszero(zz) || (valuation(zz, PP[i]) >= VD[i])
            push!(z, F(0))
          else
            push!(z, F(1))
          end
        end
        push!(S, V(z)::elem_type(V))
      end
      #S = [ V(elem_type(F)[valuation(R(EabstoE(evaluate(f(k(u)))) - 1), PP[i]) >= VD[i] ? F(0) : F(1) for i in 1:length(PP)]) for u in gens(K)]
      _T, _ = sub(V, S)
      W, w = quo(V, _T)
      if dim(W) == 0
        PP = ideal_type(R)[]
      end
    end
  end

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

  if isempty(PP)
    S = GrpAbFinGenElem[]
    Q, q = quo(Q0, S)::Tuple{GrpAbFinGen, GrpAbFinGenMap}
    Work = isdefinite(L) ? typeof(P0)[ P0 ] : typeof(P0)[]
    p = 2
    while order(Q) > 1
      while isempty(Work)
        p = next_prime(p)
        Work = ideal_type(R)[ QQ for QQ in support(p * R) if issplit(QQ) && valuation(bad, QQ) == 0 ]
      end
      P = popfirst!(Work)
      c = (q00\P)::GrpAbFinGenElem
      o = order(q(c))::fmpz
      if !isone(o)
        push!(S, c)
        Q, q = quo(Q0, S)::Tuple{GrpAbFinGen, GrpAbFinGenMap}
        push!(Gens, (P, o))
      end
    end
  else
    ll = Int(order(Q0))
    cocycle = Matrix{elem_type(W)}(undef, ll, ll)
    for i in 1:ll
      for j in 1:ll
        cocycle = zero(W)
      end
    end
    C = collect(Q0)
    ideals = ideal_type(Rabs)[ q00(C[i]) for i in 1:length(C) ]
    for i in 1:ll
      for j in 1:ll
        ij = findfirst(isequal(C[i] * C[j]), C)
        I = ideals[i] * ideals[j] * inv(ideals[ij])
        J = I * inv(a(I))
        ok, x = isprincipal(J)
        u = f(nnorm(-(ff\FacElem(nf(RR)(norm(x))))))
        x = x * u
        @assert norm(x) == 1
        y = w(V(Int[ valuation(x - 1, PP[i]) >= VD[i] ? F(0) : F(1) for i in 1:length(PP)]))
        cocycle[i, j] = y
        cocycle[j, i] = y
      end
    end

    S = Tuple{elem_type(Q0), Generic.QuotientModuleElem{elem_type(F)}}[(id(Q0), zero(W))]
    Work = isdefinite(L) ? typeof(P0)[ P0 ] : typeof(P0)[]
    p = 2
    while length(S) != order(Q0) * length(W)
      while isempty(Work)
        p = next_prime(p)
        Work = ideal_type(R)[ QQ for QQ in support(p * R) if issplit(QQ) && valuation(bad, QQ) == 0 ]
      end
      P = popfirst!(Work)
      c = q00\P
      i = findfirst(isequal(c), C)
      I = P * inv(ideals[i])
      J = I * inv(a(I))
      ok, x = isprincipal(J)
      u = f(nnorm(-(ff\FacElem(nf(RR)(norm(x))))))
      x = x * u
      @assert norm(x) == 1
      y = V(Int[ valuation(x - 1, PP[i]) >= VD[i] ? F(0) : F(1) for i in 1:length(PP)])
      idx = findfirst(isequal(P), PP)
      if idx !== nothing
        y = V(elem_type(F)[i == idx ? y[i] : y[i] + 1] for i in 1:dim(V)) #w(V([y[idx] = y[idx] + 1
      end
      elt = (c, w(y))
      elt1 = elt
      o = 1
      siz = length(S)
      while !(elt1 in S)
        j = findfirst(isequal(elt1[1]), C)
        @assert !isnothing(j)
        for l in 1:siz
          elt2 = S[l]
          k = findfirst(isequal(elt2[1]), C)
          @assert !isnothing(k)
          prod = (elt1[1] * elt2[1], elt1[2] + elt2[2] + cycycle[j, k])
          if !(prod in S)
            push!(S, prod)
          end
        end
        elt1 = (elt[1] * elt1[1], elt[2] + elt[1] + cocycle[i, j])
        o = o + 1
      end
      @assert length(S) == siz * o
      if o != 1
        push!(Gens, (P, o))
      end
    end
  end

  if isdefinite(L)
    return Gens, true, P0
  else
    return Gens, false, P0
  end
end

function representatives(G::GenusHerm)
  return genus_representatives(representative(G))
end

function _all_row_span(M)
  F = base_ring(M)
  rows = Vector{Vector{elem_type(F)}}(undef, nrows(M))
  for i in 1:nrows(M)
    rows[i] = elem_type(F)[M[i, j] for j in 1:ncols(M)]
  end
  n = nrows(M)
  it = Iterators.product([F for i in 1:n]...)
  res = Vector{Vector{elem_type(F)}}()
  for c in it
    z = Ref(c[1]) .* rows[1]
    for i in 2:n
      z = z .+ (Ref(c[i]) .* rows[i])
    end
    push!(res, z)
  end
  return res
end


text:
# These types all shadow libogg datatypes, hence they are all immutable
struct OggSyncState
    # Pointer to buffered stream data
    data::Ptr{UInt8}
    # Current allocated size of the stream buffer held in *data
    storage::Cint
    # The number of valid bytes currently held in *data; functions as the buffer head pointer
    fill::Cint
    # The number of bytes at the head of *data that have already been returned as pages;
    # functions as the buffer tail pointer
    returned::Cint

    # Synchronization state flag; nonzero if sync has not yet been attained or has been lost
    unsynced::Cint
    # If synced, the number of bytes used by the synced page's header
    headerbytes::Cint
    # If synced, the number of bytes used by the synced page's body
    bodybytes::Cint

    # zero-constructor
    OggSyncState() = new(C_NULL, 0, 0, 0, 0, 0, 0)
end

function ogg_sync_destroy(sync::OggSyncState)
    ccall((:ogg_sync_destroy,libogg), Cint, (Ref{OggSyncState},),sync)
end

struct OggPage
    # Pointer to the page header for this page
    header::Ptr{UInt8}
    # Length of the page header in bytes
    header_len::Clong
    # Pointer to the data for this page
    body::Ptr{UInt8}
    # Length of the body data in bytes
    body_len::Clong

    # zero-constructor
    OggPage() = new(C_NULL, 0, C_NULL, 0)
end

function read(page::OggPage)
    header_ptr = unsafe_wrap(Vector{UInt8}, page.header, page.header_len; own=false)
    body_ptr = unsafe_wrap(Vector{UInt8}, page.body, page.body_len; own=false)
    return vcat(header_ptr, body_ptr)
end

function show(io::IO, x::OggPage)
    write(io, "OggPage ID: $(ogg_page_serialno(x)), length $(x.body_len) bytes")
end

# This const here so that we don't use ... syntax in new()
const oss_zero_header = tuple(zeros(UInt8, 282)...)
mutable struct OggStreamState
    # Pointer to data from packet bodies
    body_data::Ptr{UInt8}
    # Storage allocated for bodies in bytes (filled or unfilled)
    body_storage::Clong
    # Amount of storage filled with stored packet bodies
    body_fill::Clong
    # Number of elements returned from storage
    body_returned::Clong

    # String of lacing values for the packet segments within the current page
    # Each value is a byte, indicating packet segment length
    lacing_vals::Ptr{Cint}
    # Pointer to the lacing values for the packet segments within the current page
    granule_vals::Ptr{Int64}
    # Total amount of storage (in bytes) allocated for storing lacing values
    lacing_storage::Clong
    # Fill marker for the current vs. total allocated storage of lacing values for the page
    lacing_fill::Clong
    # Lacing value for current packet segment
    lacing_packet::Clong
    # Number of lacing values returned from lacing_storage
    lacing_returned::Clong

    # Temporary storage for page header during encode process, while the header is being created
    header::NTuple{282,UInt8}
    # Fill marker for header storage allocation. Used during the header creation process
    header_fill::Cint

    # Marker set when the last packet of the logical bitstream has been buffered
    e_o_s::Cint
    # Marker set after we have written the first page in the logical bitstream
    b_o_s::Cint
    # Serial number of this logical bitstream
    serialno::Clong
    # Number of the current page within the stream
    pageno::Clong
    # Number of the current packet
    packetno::Int64
    # Exact position of decoding/encoding process
    granulepos::Int64

    # zero-constructor
    OggStreamState() = new(0,0,0,0,C_NULL,C_NULL,0,0,0,0,oss_zero_header,0,0,0,0,0,0,0)
end

function ogg_stream_destroy(stream::OggStreamState)
    ccall((:ogg_stream_destroy,libogg), Cint, (Ref{OggStreamState},), stream)
end

struct OggPacket
    # Pointer to the packet's data. This is treated as an opaque type by the ogg layer
    packet::Ptr{UInt8}
    # Indicates the size of the packet data in bytes. Packets can be of arbitrary size
    bytes::Clong
    # Flag indicating whether this packet begins a logical bitstream
    # 1 indicates this is the first packet, 0 indicates any other position in the stream
    b_o_s::Clong
    # Flag indicating whether this packet ends a bitstream
    # 1 indicates the last packet, 0 indicates any other position in the stream
    e_o_s::Clong

    # A number indicating the position of this packet in the decoded data
    # This is the last sample, frame or other unit of information ('granule')
    # that can be completely decoded from this packet
    granulepos::Int64
    # Sequential number of this packet in the ogg bitstream
    packetno::Int64
end
# zero-constructor
OggPacket() = OggPacket(C_NULL, 0, 0, 0, 0, 0)

function show(io::IO, x::OggPacket)
    write(io,"OggPacket ID: $(x.packetno), length $(x.bytes) bytes")
end

ext: jl


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct CylindricalTank <: Tank
        name::String
        available::Bool
        junction::Junction
        diameter::Float64
        level::Union{Nothing,Float64}
        level_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        internal::InfrastructureSystemsInternal
    end

A cylindrical water tank. Currently the only supproted tank type.

# Arguments
- `name::String`
- `available::Bool`
- `junction::Junction`
- `diameter::Float64`: constant diameter in m
- `level::Union{Nothing,Float64}`: water level in m
- `level_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `internal::InfrastructureSystemsInternal`: internal reference, do not modify
"""
mutable struct CylindricalTank <: Tank
    name::String
    available::Bool
    junction::Junction
    "constant diameter in m"
    diameter::Float64
    "water level in m"
    level::Union{Nothing,Float64}
    level_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function CylindricalTank(name, available, junction, diameter, level, level_limits, )
    CylindricalTank(name, available, junction, diameter, level, level_limits, InfrastructureSystemsInternal(), )
end

function CylindricalTank(; name, available, junction, diameter, level, level_limits, )
    CylindricalTank(name, available, junction, diameter, level, level_limits, )
end

# Constructor for demo purposes; non-functional.
function CylindricalTank(::Nothing)
    CylindricalTank(;
        name="init",
        available=true,
        junction=Junction(nothing),
        diameter=0,
        level=nothing,
        level_limits=(min=0.0, max=1000.0),
    )
end

"""Get CylindricalTank name."""
get_name(value::CylindricalTank) = value.name
"""Get CylindricalTank available."""
get_available(value::CylindricalTank) = value.available
"""Get CylindricalTank junction."""
get_junction(value::CylindricalTank) = value.junction
"""Get CylindricalTank diameter."""
get_diameter(value::CylindricalTank) = value.diameter
"""Get CylindricalTank level."""
get_level(value::CylindricalTank) = value.level
"""Get CylindricalTank level_limits."""
get_level_limits(value::CylindricalTank) = value.level_limits
"""Get CylindricalTank internal."""
get_internal(value::CylindricalTank) = value.internal


text:
using Documenter, EquationsOfState

DocMeta.setdocmeta!(EquationsOfState, :DocTestSetup, :(using EquationsOfState, EquationOfState.Collections, Unitful); recursive=true)
makedocs(;
    modules=[EquationsOfState],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Installation" => "Installation.md",
        "Manual" => [
            "Collections" => "Collections.md",
            "Nonlinear fitting" => "NonlinearFitting.md",
            "Find volume" => "Find.md",
            "Portability" => "Portability.md",
            "Interoperability" => "Python.md",
        ],
        "FAQ" => "FAQ.md",
    ],
    repo="https://github.com/MineralsCloud/EquationsOfState.jl/blob/{commit}{path}#L{line}",
    sitename="EquationsOfState.jl",
    authors="Qi Zhang <singularitti@outlook.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/MineralsCloud/EquationsOfState.jl.git"
)


text:
# Uses new formulation from limbdark paper.
include("cel_bulirsch.jl")

"""
    s2(r,b)

Compute the integral of the linear limb-darkening case.

# Examples
```julia-repl
julia> s2(0.1,0.5)
2.067294367278038
```
"""
function s2(r::T,b::T) where {T <: Real}
s_2,Eofk,Em1mKdm = s2_ell(r,b)
return s_2
end

"""
    s2_ell(r,b)

Compute the integral of the linear limb-darkening case, s_2 from
starry (which is S_1 in Agol & Luger) and also return two elliptic
integrals, E(k) and (E(m)-(1-m)K(m))/m, where m = k^2.

# Examples
```julia-repl
julia> s2_ell(0.1,0.5)
(2.067294367278038, 1.4726447391114554, 0.8111640472029077)
```
"""
function s2_ell(r::T,b::T) where {T <: Real}
if T == Float64
  third = third_float
elseif T == BigFloat
  third = third_big
else
  third = convert(T,1)/convert(T,3)
end
# For now, just compute linear component:
Eofk = zero(T)
Em1mKdm = zero(T)
Lambda1 = zero(T)
if b >= 1.0+r ||  r == 0.0
  # No occultation:
  Lambda1 = zero(T)  # Case 1
elseif b <= r-1.0
  # Full occultation:
  Lambda1 = zero(T)  # Case 11
else
  if b == 0
    Lambda1 = -2pi*sqrt(1.0-r^2)^3 # Case 10
    Eofk = .5*pi
    Em1mKdm = .25*pi
  elseif b==r
    if r == 0.5
#      Lambda1 = pi-4/3 - 2*(b-0.5) +6*(r-0.5) # Case 6; I've added in analytic first derivaties.
      Lambda1 = pi-4*third - 2*(b-0.5) +6*(r-0.5) # Case 6; I've added in analytic first derivaties.
      Eofk = one(T)
      Em1mKdm = one(T)
    elseif r < 0.5
      m = 4r^2
      Eofk = cel_bulirsch(m,one(T),one(T),1-m)
      Em1mKdm = cel_bulirsch(m,one(T),one(T),zero(T))
#      Lambda1 = pi+2/3*((2m-3)*Eofk - m*Em1mKdm)+ # Case 5
      Lambda1 = pi+2*third*((2m-3)*Eofk - m*Em1mKdm)+ # Case 5
        (b-r)*4*r*(Eofk-2*Em1mKdm)  # Adding in first derivative
    else
      m = 4r^2; minv = inv(m)
      Eofk = cel_bulirsch(minv,one(T),one(T),1-minv)
      Em1mKdm = cel_bulirsch(minv,one(T),one(T),zero(T))
#      Lambda1 = pi+1/(3*r)*(-m*Eofk + (2m-3)*Em1mKdm) - # Case 7
      Lambda1 = pi+third/r*(-m*Eofk + (2m-3)*Em1mKdm) - # Case 7
        (b-r)*2*(2*Eofk-Em1mKdm) # Adding in first derivative
    end
  else
    onembpr2 = (1-b-r)*(1+b+r); onembmr2=(r-b+1)*(1-r+b); fourbr = 4b*r; fourbrinv = inv(fourbr)
    k2 = onembpr2*fourbrinv+1
    if (b+r) > 1.0 # k^2 < 1, Case 2, Case 8
      kc2 = -onembpr2*fourbrinv; kc = sqrt(kc2)
      Piofk,Eofk,Em1mKdm = cel_bulirsch(k2,kc,(b-r)^2*kc2,zero(T),one(T),one(T),3*kc2*(b-r)*(b+r),kc2,zero(T))
#      Lambda1 = onembmr2*(Piofk+ (-3+6r^2+2*b*r)*Em1mKdm-fourbr*Eofk)/(3*sqrt(b*r))
      Lambda1 = onembmr2*(Piofk+ (-3+6r^2+2*b*r)*Em1mKdm-fourbr*Eofk)*third/sqrt(b*r)
    elseif (b+r) < 1.0  # k^2 > 1, Case 3, Case 9
      onembmr2inv=inv(onembmr2); k2inv = inv(k2); kc2 =onembpr2*onembmr2inv; kc = sqrt(kc2)
      bmrdbpr = (b-r)/(b+r)
      mu = 3bmrdbpr*onembmr2inv
      p = bmrdbpr^2*onembpr2*onembmr2inv
      Piofk,Eofk,Em1mKdm = cel_bulirsch(k2inv,kc,p,1+mu,one(T),one(T),p+mu,kc2,zero(T))
#      Lambda1 = 2*sqrt(onembmr2)*(onembpr2*Piofk -(4-7r^2-b^2)*Eofk)/3
      Lambda1 = 2*sqrt(onembmr2)*(onembpr2*Piofk -(4-7r^2-b^2)*Eofk)*third
    else
      # b+r = 1 or k^2=1, Case 4 (extending r up to 1)
#      Lambda1 = 2*acos(1.0-2.0*r)-2pi*convert(T,r>.5) -(4/3*(3+2r-8r^2)+
      Lambda1 = 2*acos(1.0-2.0*r)-2pi*convert(T,r>.5) -(4*third*(3+2r-8r^2)+
          8*(r+b-1)*r)*sqrt(r*(1-r)) # Adding in first derivatives
      Eofk = one(T)
      Em1mKdm = one(T)
    end
  end
end
#flux = ((1.0-convert(T,r>b))*2pi-Lambda1)/3
flux = ((1.0-convert(T,r>b))*2pi-Lambda1)*third
return flux,Eofk,Em1mKdm
end

"""
    s2!(r,b,s2_grad)

Computes the linear limb-darkening case, as well as the gradient,
s2_grad=[ds_2/dr,ds_2/db] is a pre-allocated two-element array.
Returns s2 and complete elliptic integrals, E(k) and
(E(m)-(1-m)K(m))/m, where m = k^2.

# Example
```julia-repl
julia> s2_grad=[0.0,0.0]
julia> s2!(0.1,0.5,s2_grad)
(2.067294367278038, 1.4726447391114554, 0.8111640472029077)
julia> s2_grad
2-element Array{Float64,1}:
 -0.53988
  0.0182916
```
"""
function s2!(r::T,b::T,s2_grad::Array{T,1}) where {T <: Real}
if T == Float64
  third = third_float
elseif T == BigFloat
  third = third_big
else
  third = convert(T,1)/convert(T,3)
end
Eofk = zero(T)
Em1mKdm = zero(T)
Lambda1 = zero(T)
fill!(s2_grad,zero(T))
if b >= 1.0+r ||  r == 0.0
  # No occultation:
  Lambda1 = zero(T)  # Case 1
elseif b <= r-1.0
  # Full occultation:
  Lambda1 = zero(T)  # Case 11
else
  if b == 0
    sqrt1mr2 = sqrt(1.0-r^2)
    Lambda1 = -2pi*sqrt1mr2^3 # Case 10
    s2_grad[1] = -2pi*r*sqrt1mr2 # dLambda/dr (dLambda/db= 0)
    Eofk = .5*pi
    Em1mKdm = .25*pi
  elseif b==r
    if r == 0.5
#      Lambda1 = pi-4/3 # Case 6; added in analytic first derivaties.
      Lambda1 = pi-4*third # Case 6; added in analytic first derivaties.
      s2_grad[1] = -2      # dLambda/dr
#      s2_grad[2] =  2/3  # dLambda/db
      s2_grad[2] =  2*third  # dLambda/db
      Eofk = one(T)
      Em1mKdm = one(T)
    elseif r < 0.5
      m = 4r^2
      Eofk = cel_bulirsch(m,one(T),one(T),1-m)
      Em1mKdm = cel_bulirsch(m,one(T),one(T),zero(T))
#      Lambda1 = pi+2/3*((2m-3)*Eofk - m*Em1mKdm) # Case 5
      Lambda1 = pi+2*third*((2m-3)*Eofk - m*Em1mKdm) # Case 5
      s2_grad[1] = -4*r*Eofk      # Adding in first derivative dLambda/dr
#      s2_grad[2] = -4*r/3*(Eofk-2*Em1mKdm) # Adding in first derivative dLambda/db
      s2_grad[2] = -4*r*third*(Eofk-2*Em1mKdm) # Adding in first derivative dLambda/db
    else
      m = 4r^2; minv = inv(m); kc = sqrt(1.0-minv)
      Eofk = cel_bulirsch(minv,one(T),one(T),1-minv)
      Em1mKdm = cel_bulirsch(minv,one(T),one(T),zero(T))
#      Lambda1 = pi+1/(3*r)*(-m*Eofk + (2m-3)*Em1mKdm)  # Case 7
      Lambda1 = pi+third/r*(-m*Eofk + (2m-3)*Em1mKdm)  # Case 7
      s2_grad[1] = -2*Em1mKdm # dLambda/dr
#      s2_grad[2] =  2/3*(2*Eofk-Em1mKdm)  # dLambda/db
      s2_grad[2] =  2*third*(2*Eofk-Em1mKdm)  # dLambda/db
    end
  else
    onembpr2 = (1-r-b)*(1+r+b); onembmr2=(r+1-b)*(1-r+b); fourbr = 4b*r; fourbrinv=inv(fourbr)
    k2 = onembpr2*fourbrinv+1
    if (b+r) > 1.0 # k^2 < 1, Case 2, Case 8
      kc2 = -onembpr2*fourbrinv; kc = sqrt(kc2); sqbr=sqrt(b*r); sqbrinv = inv(sqbr)
      Piofk,Eofk,Em1mKdm = cel_bulirsch(k2,kc,(b-r)^2*kc2,zero(T),one(T),one(T),3*kc2*(b-r)*(b+r),kc2,zero(T))
#      Lambda1 = onembmr2*(Piofk+ (-3+6r^2+2*b*r)*Em1mKdm-fourbr*Eofk)*sqbrinv/3
      Lambda1 = onembmr2*(Piofk+ (-3+6r^2+2*b*r)*Em1mKdm-fourbr*Eofk)*sqbrinv*third
      s2_grad[1] = -2r*onembmr2*Em1mKdm*sqbrinv
#      s2_grad[2] = 2r*onembmr2*(-Em1mKdm+2*Eofk)*sqbrinv/3
      s2_grad[2] = 2r*onembmr2*(-Em1mKdm+2*Eofk)*sqbrinv*third
    elseif (b+r) < 1.0  # k^2 > 1, Case 3, Case 9
      onembmr2inv = inv(onembmr2); k2inv = inv(k2); kc2 =onembpr2*onembmr2inv; kc = sqrt(kc2)
      bmrdbpr = (b-r)/(b+r);
      mu = 3bmrdbpr*onembmr2inv
      p = bmrdbpr^2*onembpr2*onembmr2inv
#      println("calling cel with: ",k2inv,kc,p,1+mu,one(T),one(T),p+mu,kc2,zero(T))
      Piofk,Eofk,Em1mKdm = cel_bulirsch(k2inv,kc,p,1+mu,one(T),one(T),p+mu,kc2,zero(T))
      sqonembmr2 = sqrt(onembmr2)
#      Lambda1 = 2*sqonembmr2*(onembpr2*Piofk -(4-7r^2-b^2)*Eofk)/3
      Lambda1 = 2*sqonembmr2*(onembpr2*Piofk -(4-7r^2-b^2)*Eofk)*third
      s2_grad[1] = -4*r*sqonembmr2*Eofk
#      s2_grad[2] = -4*r/3*sqonembmr2*(Eofk - 2*Em1mKdm)
      s2_grad[2] = -4*r*third*sqonembmr2*(Eofk - 2*Em1mKdm)
    else
      # b+r = 1 or k^2=1, Case 4 (extending r up to 1)
#      Lambda1 = 2*acos(1.0-2.0*r)-4/3*(3+2r-8r^2)*sqrt(r*(1-r))-2pi*convert(T,r>0.5)
      Lambda1 = 2*acos(1.0-2.0*r)-4*third*(3+2r-8r^2)*sqrt(r*(1-r))-2pi*convert(T,r>0.5)
      s2_grad[1] = -8*r*sqrt(r*(1-r))
#      s2_grad[2] = -s2_grad[1]/3
      s2_grad[2] = -s2_grad[1]*third
      Eofk = one(T)
      Em1mKdm = one(T)
    end
  end
end
#flux = ((1.0-convert(T,r>b))*2pi-Lambda1)/3
flux = ((1.0-convert(T,r>b))*2pi-Lambda1)*third
return flux,Eofk,Em1mKdm
end

#"""
#    s2!(trans)
#
#Computes the linear limb-darkening case, as well as the gradient,
#trans.s2_grad=[ds_2/dr,ds_2/db] is a pre-allocated two-element array.
#Returns trans.s2, and computes the complete elliptic integrals,
#trans.Eofk, E(k), and trans.Em1mKdm, (E(m)-(1-m)K(m))/m, where m = k^2.
#
#NOTE:  This is only meant to be run in tandem with `transit_poly_struct.jl`,
#which computes the necessary elements of the structure first.
#
## Example
#```julia-repl
#julia> trans=transit_init(0.1,0.5,[0.5,0.5],true)
#julia> s2!(trans)
#julia> trans.sn[2]
#julia> trans.s2_grad
#2-element Array{Float64,1}:
# -0.53988
#  0.0182916
#```
#"""
#function s2!(t::Transit_Struct{T}) where {T <: Real}
#r=t.r; b=t.b
## Elliptic integrals:
#t.Eofk = zero(T)
#t.Em1mKdm = zero(T)
## Lambda1 = 3\pi \Lambda_1 from starry/limbdark papers:
#Lambda1 = zero(T)
#if t.grad
#  fill!(t.s2_grad,zero(T))
#end
#if b >= 1.0+r ||  r == 0.0
#  # No occultation:
#  Lambda1 = zero(T)  # Case 1
#elseif b <= r-1.0
#  # Full occultation:
#  Lambda1 = zero(T)  # Case 11
#else
#  if b == 0
#    Lambda1 = -2pi*t.sqrt1mr2^3 # Case 10
#    if t.grad
#      t.s2_grad[1] = -2pi*r*t.sqrt1mr2 # dLambda/dr (dLambda/db= 0)
#    end
#    t.Eofk = 0.5*pi
#    t.Em1mKdm = 0.25*pi
#  elseif b==r
#    if r == 0.5
#      Lambda1 = pi-4*t.third # Case 6; added in analytic first derivaties.
#      if t.grad
#        t.s2_grad[1] = -2      # dLambda/dr
#        t.s2_grad[2] =  t.twothird  # dLambda/db
#      end
#      t.Eofk = one(T)
#      t.Em1mKdm = one(T)
#    elseif r < 0.5
#      m = 4r^2
#      t.Eofk = cel_bulirsch(m,one(T),one(T),1-m)
#      t.Em1mKdm = cel_bulirsch(m,one(T),one(T),zero(T))
#      Lambda1 = pi+2*t.third*((2m-3)*t.Eofk - m*t.Em1mKdm) # Case 5
#      if t.grad
#        t.s2_grad[1] = -4*r*t.Eofk      # Adding in first derivative dLambda/dr
#        t.s2_grad[2] = -4*r*t.third*(t.Eofk-2*t.Em1mKdm) # Adding in first derivative dLambda/db
#      end
#    else
#      m = 4r^2; minv = inv(m)
#      t.Eofk = cel_bulirsch(minv,one(T),one(T),1-minv)
#      t.Em1mKdm = cel_bulirsch(minv,one(T),one(T),zero(T))
#      Lambda1 = pi+t.third/r*(-m*t.Eofk + (2m-3)*t.Em1mKdm)  # Case 7
#      if t.grad
#        t.s2_grad[1] = -2*t.Em1mKdm # dLambda/dr
#        t.s2_grad[2] =  2*t.third*(2*t.Eofk-t.Em1mKdm)  # dLambda/db
#      end
#    end
#  else
#    t.onembpr2 = (1-b-r)*(1+b+r); t.onembmr2=(r-b+1)*(1-r+b); t.fourbr = 4b*r; t.fourbrinv = inv(t.fourbr)
#    t.k2 = t.onembpr2*t.fourbrinv+1
#    if (b+r) > 1.0 # k^2 < 1, Case 2, Case 8
#      t.kc2 = -t.onembpr2*t.fourbrinv; t.kc = sqrt(t.kc2); t.sqbr=sqrt(b*r); t.sqbrinv = inv(t.sqbr)
#      Piofk,t.Eofk,t.Em1mKdm = cel_bulirsch(t.k2,t.kc,(b-r)^2*t.kc2,zero(T),one(T),one(T),3*t.kc2*(b-r)*(b+r),t.kc2,zero(T))
#      Lambda1 = t.onembmr2*(Piofk+ (-3+6r^2+2*b*r)*t.Em1mKdm-t.fourbr*t.Eofk)*t.sqbrinv*t.third
#      if t.grad
#        t.s2_grad[1] = -2r*t.onembmr2*t.Em1mKdm*t.sqbrinv
#        t.s2_grad[2] = 2r*t.onembmr2*(-t.Em1mKdm+2*t.Eofk)*t.sqbrinv*t.third
#      end
#    elseif (b+r) < 1.0  # k^2 > 1, Case 3, Case 9
#      t.onembmr2inv=inv(t.onembmr2); t.k2inv = inv(t.k2); t.kc2 = t.onembpr2*t.onembmr2inv; t.kc = sqrt(t.kc2)
#      bmrdbpr = (b-r)/(b+r);
#      mu = 3bmrdbpr*t.onembmr2inv
#      p = bmrdbpr^2*t.onembpr2*t.onembmr2inv
#      t.k2inv = inv(t.k2)
##      println("calling cel with: ",t.k2inv,t.kc,p,1+mu,one(T),one(T),p+mu,t.kc2,zero(T))
#      Piofk,t.Eofk,t.Em1mKdm = cel_bulirsch(t.k2inv,t.kc,p,1+mu,one(T),one(T),p+mu,t.kc2,zero(T))
#      t.sqonembmr2 = sqrt(t.onembmr2)
#      Lambda1 = 2*t.sqonembmr2*(t.onembpr2*Piofk -(4-7r^2-b^2)*t.Eofk)*t.third
#      if t.grad
#        t.s2_grad[1] = -4*r*t.sqonembmr2*t.Eofk
#        t.s2_grad[2] = -4*r*t.third*t.sqonembmr2*(t.Eofk - 2*t.Em1mKdm)
#      end
#    else
#      # b+r = 1 or k^2=1, Case 4 (extending r up to 1)
#      t.sqr1mr = sqrt(r*(1-r))
#      Lambda1 = 2*acos(1.0-2.0*r)-4*t.third*(3+2r-8r^2)*t.sqr1mr-2pi*convert(T,r>0.5)
#      if t.grad
#        t.s2_grad[1] = -8*r*t.sqr1mr
#        t.s2_grad[2] = -t.s2_grad[1]*t.third
#      end
#      t.Eofk = one(T)
#      t.Em1mKdm = one(T)
#    end
#  end
#end
##flux = ((1.0-convert(T,r>b))*2pi-Lambda1)/3
#t.sn[2] = ((1.0-convert(T,r>b))*2pi-Lambda1)*t.third
#return
#end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
######################################################################
# general interface
######################################################################

export
    # general interface
    AbstractBasis,
    domain,
    collocation_points,
    degf,
    basis_matrix,
    InterpolatedFunction,
    collocation_values,
    # bases
    ChebyshevBasis,
    IntervalAB,
    LinearInterpolation,
    ExtrapolateLevel

abstract AbstractBasis

"""
Return the domain of the basis.
"""
function domain end

"A range of `n` linearly spaced elements in the domain of `basis`."
function linspace(basis::AbstractBasis, len)
    dom = domain(basis)
    linspace(dom.lo, dom.hi, len)
end

"""
Return the collocation points (if any) for a basis.
"""
function collocation_points end

"""
Degress of freedom (number of coefficients).
"""
degf(basis::AbstractBasis) = length(collocation_points(basis))

"""
Return the basis matrix, ie the functions in the basis evaluated at
the given points, one row for each point. Points default to
`collocation_points`.
"""
function basis_matrix end

"""
Interpolated function using `basis`, with coefficients `α`.
"""
immutable InterpolatedFunction{S, T}
    basis::S
    α::T
end

*(basis::AbstractBasis, α) = InterpolatedFunction(basis, α)

@forward InterpolatedFunction.basis domain, collocation_points, degf

function show(io::IO, ipf::InterpolatedFunction)
    print(io, "Interpolated function on ", ipf.basis)
end

"""
Evaluate the approximation on `basis` with coefficients `α` at at `x`.

Notes:
1. methods should be defined for various bases.
2. the callable interface is recommended for use outside this module.
"""
function _evaluate(basis, α, x)
    first(basis_matrix(basis, [x]) * α)
end

(f::InterpolatedFunction)(x) = _evaluate(f.basis, f.α, x)

"""
Evaluate the interpolated function at the collocation points.
"""
function collocation_values(ipf::InterpolatedFunction)
    ipf.(collocation_points(ipf.basis))
end

"""
Interpolate points in `y`, assumed to be values at the collocation
points, using `basis`.
"""
function \(basis::AbstractBasis, y::AbstractVector)
    # NOTE: this is just a fallback method, not necessarily efficient.
    @assert degf(basis) == length(y)
    InterpolatedFunction(basis, basis_matrix(basis) \ y)
end

"""
Interpolate `f`, evaluated at the collocation points, using `basis`.
"""
function \(basis::AbstractBasis, f)
    basis \ f.(collocation_points(basis))
end

zeros(basis::AbstractBasis) = basis \ zeros(degf(basis))

ones(basis::AbstractBasis) = basis \ ones(degf(basis))

# plot recipe for displaying interpolated functions
@recipe function f(ipf::InterpolatedFunction, label = "function", N = 100)
    @unpack basis, α = ipf
    dom = domain(basis)
    z = collect(linspace(dom.lo, dom.hi, N))
    @series begin
        seriestype := :path
        label := label
        z, ipf.(z)
    end
    @series begin
        primary := false
        seriestype := :scatter
        collocation_points(ipf), collocation_values(ipf)
    end
end

######################################################################
# Chebyshev basis on [-1,1]
######################################################################

immutable ChebyshevBasis <: AbstractBasis
    n::Int
    stretch::Bool
    extrapolate::Bool
    function ChebyshevBasis(n; stretch=false, extrapolate=false)
        new(n, stretch, extrapolate)
    end
end

size(b::ChebyshevBasis) = (b.n,)

domain(b::ChebyshevBasis, T=Float64) = -one(T)..one(T)

function collocation_points(b::ChebyshevBasis)
    zs = [cos((i-0.5)*pi/b.n) for i in b.n:-1:1]
    if b.stretch
        scale = 1/zs[end]
        zs .= zs .* scale
    end
    zs
end

function basis_matrix{T}(b::ChebyshevBasis,
                         x::AbstractVector{T}=collocation_points(b))
    if !b.extrapolate
        @assert all(abs(x) ≤ one(T) + √eps(T) for x in x)
    end
    A = Array{T}(length(x), b.n)
    for i in 1:b.n
        A[:,i] .= i == 1 ? one(T) : i == 2 ? x : (2.*x.*A[:, i-1]-A[:, i-2])
    end
    A
end

######################################################################
# IntervalAB
######################################################################

"""
Transform an inner basis to an interval.

Transforming to the inner basis is with `x ↦ fma(scale, shift)`.
"""
immutable IntervalAB{T,S <: AbstractBasis} <: AbstractBasis
    domain::Interval{T}
    inner_basis::S
    scale::T
    shift::T
    function IntervalAB(dom, inner_basis)
        @assert isfinite(dom)
        inner_dom = domain(inner_basis)
        scale = diam(inner_dom)/diam(dom)
        new(dom, inner_basis, scale, mid(inner_dom) - scale*mid(dom))
    end
end

function IntervalAB{T,S}(domain::Interval{T}, inner_basis::S)
    IntervalAB{T,S}(domain, inner_basis)
end

_map_to_inner(b::IntervalAB, x) = fma.(x, b.scale, b.shift)

_map_from_inner(b::IntervalAB, x) = fma.(x, 1/b.scale, -b.shift/b.scale)

domain(b::IntervalAB) = b.domain

function collocation_points(b::IntervalAB)
    _map_from_inner(b, collocation_points(b.inner_basis))
end

basis_matrix(b::IntervalAB) = basis_matrix(b.inner_basis)

basis_matrix(b::IntervalAB, x) = basis_matrix(b.inner_basis, _map_to_inner(b, x))

function show(io::IO, b::IntervalAB)
    print(io, b.inner_basis)
    print(io, " on ")
    print(io, b.domain)
end

######################################################################
# linear interpolation
######################################################################

immutable LinearInterpolation{T <: AbstractVector} <: AbstractBasis
    nodes::T
    function LinearInterpolation(nodes)
        @assert issorted(nodes) "Nodes need to be sorted"
        new(nodes)
    end
end

LinearInterpolation{T}(nodes::T) = LinearInterpolation{T}(nodes)

collocation_points(li::LinearInterpolation) = li.nodes

domain(li::LinearInterpolation) = li.nodes[1]..li.nodes[end]

basis_matrix(li::LinearInterpolation) = speye(length(li.nodes))

function \(li::LinearInterpolation, y::AbstractVector)
    @assert degf(li) == length(y)
    InterpolatedFunction(li, y)
end

function _evaluate(li::LinearInterpolation, α, x)
    @unpack nodes = li
    (nodes[1] ≤ x ≤ nodes[end]) || error(DomainError())
    index = searchsortedlast(nodes, x)
    @assert index ≠ 0 "internal error"
    # println("index = $index, node = $(nodes[index])")
    if x == nodes[index]
        α[index]
    else
        p = (x - nodes[index]) / (nodes[index+1] - nodes[index])
        α[index]*(1-p) + α[index+1]*p
    end
end

function collocation_values{S <: LinearInterpolation, T}(ipf::InterpolatedFunction{S,T})
    ipf.α
end

######################################################################
# ExtrapolateLevel
######################################################################

immutable ExtrapolateLevel{T} <: AbstractBasis
    inner_basis::T
end

@forward ExtrapolateLevel.inner_basis domain, collocation_points, basis_matrix, degf

function _evaluate(basis::ExtrapolateLevel, α, x)
    @unpack inner_basis = basis
    dom = domain(inner_basis)
    _evaluate(inner_basis, α, clamp(x, dom.lo, dom.hi))
end


text:
module SyntheticImplicitFeedback

export Feature, Rule, accumulate, generate

type Feature
    name::UTF8String
    default::Any
    random::Function
end

type Rule
    # return bool
    f::Function

    # if f returns true, accumulative CTR is lifted p%
    p::Float64
end

function accumulate(sample::Dict, rules::Array{Rule,1})
    ctr = 0.0
    for rule in rules
        if rule.f(sample)
            ctr += rule.p
        end
    end
    ctr
end

function generate(sample::Dict, rules::Array{Rule,1})
    rand() <= accumulate(sample, rules)
end

function generate(samples::Array{Dict,1}, rules::Array{Rule,1})
    map(samples) do sample
        generate(sample, rules)
    end
end

end # module


text:
using LabelledArrays

export chaoticmaps

# """
#     eom_linearmap1(x, p, n) → SVector{2}

# # References
# Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended
# Granger causality." Physics Letters A 324.1 (2004): 26-35
# """
# function eom_linearmap1(x, p, t)
#     c = p[1]
#     x, y = (x...,)
#     dx = 3.4*x*(t - 1)* (1 - x^2*(t - 1)) * exp((-x^2)*(t - 1)) + 0.8*x*(t - 2)
#     dy = 3.4*y*(t - 1)* (1 - y^2*(t - 1)) * exp((-y^2)*(t - 1)) + 0.5*y*(t - 2) + c*x*(t - 2)
#     return SVector{2}(dx, dy)
# end

# """
#     linearmap1(u₀, c) → DiscreteDynamicalSystem

# # References
# Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended Granger causality." Physics Letters A 324.1 (2004): 26-35
# """
# function linearmap1(u₀, c)
#     p = @LArray [c] (:c)
#     logistic_system = DiscreteDynamicalSystem(eom_linearmap1, u₀, p)
#     return logistic_system
# end

# """
#     linearmap1(;u₀ = [1, rand(2)], c = 2.0) → DiscreteDynamicalSystem

# # References
# Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended Granger causality." Physics Letters A 324.1 (2004): 26-35
# """
# linearmap1(;u₀ = rand(2), c = 0.5) = linearmap1(u₀, c)

"""
    eom_chaoticmaps(x, p, n) → SVector{3}

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended
Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
function eom_chaoticmaps(x, p, t)
    c = p[1]
    x, y, z = (x...)
    t = t + 2
    dx = 3.4*x*(t - 1) * (1 - x^2*(t - 1)) * exp((-x^2)*(t - 1))
    dy = 3.4*y*(t - 1) * (1 - y^2*(t - 1)) * exp((-y^2)*(t - 1)) + 0.5*x*(t - 1)
    dz = 3.4*z*(t - 1) * (1 - z^2*(t - 1)) * exp((-z^2)*(t - 1)) + c*x*(t - 1) + 0.3*y*(t - 1)

    return SVector{3}(dx, dy, dz)
end

function chaoticmaps(u₀, c)
    p = [c]
    cmap = DiscreteDynamicalSystem(eom_chaoticmaps, u₀, p)
    return cmap
end

"""
    chaoticmaps(;u₀ = rand(3), c = 2.0) → DiscreteDynamicalSystem

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
chaoticmaps(;u₀ = rand(3), c = 0.5) = chaoticmaps(u₀, c)


text:
synth_paletted_imgs = [
    ("RGB_paletted", IndirectArray(rand(1:100, 127, 257), rand(RGB{N0f8}, 100))),
    ("BRG_paletted", IndirectArray(rand(1:100, 127, 257), rand(BGR{N0f8}, 100))),
    ("ARGB_paletted", IndirectArray(rand(1:100, 127, 257), rand(ARGB{N0f8}, 100))),
    ("ABGR_paletted", IndirectArray(rand(1:100, 127, 257), rand(ABGR{N0f8}, 100))),
    ("RGBA_paletted", IndirectArray(rand(1:100, 127, 257), rand(RGBA{N0f8}, 100))),
    ("RGB_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(RGB{Float64}, 100))),
    ("BRG_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(BGR{Float64}, 100))),
    ("ARGB_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(ARGB{Float64}, 100))),
    ("ABGR_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(ABGR{Float64}, 100))),
    ("RGBA_paletted_float", IndirectArray(rand(1:100, 127, 257), rand(RGBA{Float64}, 100))),
]

expected_img(x::Matrix{<:TransparentColor}) = RGBA{N0f8}.(x)
expected_img(x::Matrix{<:AbstractRGB}) = RGB{N0f8}.(x)

@testset "synthetic paletted images" begin
    for expand_paletted in (true, false)
        @testset "expand_paletted = $(expand_paletted)" begin
            for (case, image) in synth_paletted_imgs
                @testset "$(case)" begin
                    expected = expected_img(collect(image))
                    fpath = joinpath(PNG_TEST_PATH, "test_img_$(case).png")
                    @testset "write" begin
                        open(io -> PNGFiles.save(io, image), fpath, "w")
                        @test PNGFiles.save(fpath, image) === nothing
                    end
                    @testset "read" begin
                        global read_in_pngf = PNGFiles.load(fpath, expand_paletted=expand_paletted)
                        @test typeof(read_in_pngf) <: AbstractMatrix
                    end
                    @testset "compare" begin
                        @test all(expected .≈ read_in_pngf)
                    end
                    global read_in_immag = _standardize_grayness(ImageMagick.load(fpath))
                    @testset "$(case): ImageMagick read type equality" begin
                        # The lena image is Grayscale saved as RGB...
                        @test eltype(_standardize_grayness(read_in_pngf)) == eltype(read_in_immag)
                    end
                    @testset "$(case): ImageMagick read values equality" begin
                        imdiff_val = imdiff(read_in_pngf, read_in_immag)
                        onfail(@test imdiff_val < 0.01) do
                            PNGFiles._inspect_png_read(fpath)
                            _add_debugging_entry(fpath, case, imdiff_val)
                        end
                    end
                    path, ext = splitext(fpath)
                    newpath = path * "_new" * ext
                    PNGFiles.save(newpath, read_in_pngf)
                    @testset "$(case): IO is idempotent" begin
                        @test all(read_in_pngf .≈ PNGFiles.load(newpath, expand_paletted=expand_paletted))
                        @test all(read_in_pngf .≈ open(io -> PNGFiles.load(io, expand_paletted=expand_paletted), newpath))
                    end
                end
            end
        end
    end
end

@testset "palleted image with fewer transparency values than color values" begin
    expected_palette = RGBA{N0f8}[
        RGBA{N0f8}(0.0, 0.0, 0.0, 0.0),
        RGBA{N0f8}(0.282, 0.133, 0.169, 0.024),
        RGBA{N0f8}(0.502, 0.431, 0.459, 0.337),
        RGBA{N0f8}(0.49, 0.427, 0.478, 0.51),
        RGBA{N0f8}(0.282, 0.133, 0.169, 0.208),
        RGBA{N0f8}(0.0, 0.0, 0.0, 0.004),
        RGBA{N0f8}(0.714, 0.518, 0.153, 0.722),
        RGBA{N0f8}(0.6, 0.49, 0.322, 0.953),
        RGBA{N0f8}(0.298, 0.18, 0.208, 0.718),
        RGBA{N0f8}(0.196, 0.11, 0.11, 0.408),
        RGBA{N0f8}(0.0, 0.0, 0.0, 0.039),
        RGBA{N0f8}(0.655, 0.765, 0.18, 0.227),
        RGBA{N0f8}(0.592, 0.686, 0.18, 0.212),
        RGBA{N0f8}(0.961, 0.839, 0.153, 0.957),
        RGBA{N0f8}(0.58, 0.412, 0.11, 0.831),
        RGBA{N0f8}(0.0, 0.0, 0.0, 0.345),
        RGBA{N0f8}(0.0, 0.0, 0.0, 0.188),
        RGBA{N0f8}(0.0, 0.0, 0.0, 0.02),
        RGBA{N0f8}(0.667, 0.757, 0.208, 0.749),
        RGBA{N0f8}(0.616, 0.69, 0.196, 0.988),
        RGBA{N0f8}(0.584, 0.365, 0.169, 0.996),
        RGBA{N0f8}(0.427, 0.439, 0.208, 0.573),
        RGBA{N0f8}(0.537, 0.533, 0.18, 0.341),
        RGBA{N0f8}(0.51, 0.482, 0.196, 0.122),
        RGBA{N0f8}(0.545, 0.678, 0.196, 0.804),
        RGBA{N0f8}(0.208, 0.133, 0.082, 0.694),
        RGBA{N0f8}(0.298, 0.239, 0.133, 0.569),
        RGBA{N0f8}(0.082, 0.0, 0.0, 0.275),
        RGBA{N0f8}(0.0, 0.0, 0.0, 0.008),
        RGBA{N0f8}(0.518, 0.627, 0.196, 0.659),
        RGBA{N0f8}(0.427, 0.522, 0.153, 0.8),
        RGBA{N0f8}(0.412, 0.494, 0.133, 0.765),
        RGBA{N0f8}(0.239, 0.18, 0.082, 0.769),
        RGBA{N0f8}(0.196, 0.0, 0.082, 0.635),
        RGBA{N0f8}(0.0, 0.0, 0.0, 0.027),
        RGBA{N0f8}(0.082, 0.0, 0.0, 0.035),
        RGBA{N0f8}(0.196, 0.22, 0.0, 0.082),
        RGBA{N0f8}(0.082, 0.11, 0.0, 0.239),
        RGBA{N0f8}(0.208, 0.231, 0.0, 0.204),
        RGBA{N0f8}(0.133, 0.0, 0.0, 0.098),
        RGBA{N0f8}(0.133, 0.0, 0.0, 0.4),
        RGBA{N0f8}(0.0, 0.0, 0.0, 0.035),
        RGBA{N0f8}(0.796, 0.659, 0.22, 1.0), # <- Implicitly opaque
        RGBA{N0f8}(0.553, 0.698, 0.208, 1.0),
        RGBA{N0f8}(0.478, 0.561, 0.251, 1.0),
        RGBA{N0f8}(0.376, 0.314, 0.169, 1.0),
        RGBA{N0f8}(1.0, 1.0, 1.0, 1.0),
    ]
    img = PNGFiles.load("test_images/paletted_implicit_transparency.png")
    @testset "palette equality" begin
        @test img.values.parent == expected_palette
    end

    @testset begin "idempotency"
        new_path = "test_images/paletted_implicit_transparency_new.png"
        PNGFiles.save(new_path, img)
        new_img = PNGFiles.load(new_path)
        @test img == new_img
        Base.rm(new_path)
    end
end