text:
################################################################################
#
#  NfOrd/ResidueRingMultGrp.jl : Multiplicative group of Residue Rings
#
################################################################################

export multiplicative_group, multiplicative_group_generators

################################################################################
#
#  High Level Interface
#
################################################################################

@doc Markdown.doc"""
    multiplicative_group(Q::NfOrdQuoRing) -> GrpAbFinGen, Map{GrpAbFinGen, NfOrdQuoRing}
    unit_group(Q::NfOrdQuoRing) -> GrpAbFinGen, Map{GrpAbFinGen, NfOrdQuoRing}

Returns the unit group of $Q$ as an abstract group $A$ and
an isomorphism map $f \colon A \to Q^\times$.
"""
function multiplicative_group(Q::NfOrdQuoRing)
  if !isdefined(Q, :multiplicative_group)
    if is_maximal_known_and_maximal(base_ring(Q))
      G, GtoQ = _multgrp(Q)
    else
      G, GtoQ = _multgrp_non_maximal(Q)
    end
    Q.multiplicative_group = GtoQ
  end
  mQ = Q.multiplicative_group::GrpAbFinGenToNfOrdQuoRingMultMap
  return domain(mQ), mQ
end

unit_group(Q::NfOrdQuoRing) = multiplicative_group(Q)

@doc Markdown.doc"""
    multiplicative_group_generators(Q::NfOrdQuoRing) -> Vector{NfOrdQuoRingElem}

Return a set of generators for $Q^\times$.
"""
function multiplicative_group_generators(Q::NfOrdQuoRing)
  return multiplicative_group(Q).generators
end

# Factors Q.ideal, the result is saved in Q.factor
function factor(Q::NfOrdQuoRing)
  if !isdefined(Q, :factor)
    Q.factor = factor(Q.ideal)
  end
  return Q.factor
end

################################################################################
#
#  Internals
#
################################################################################

@doc Markdown.doc"""
    _multgrp(Q::NfOrdQuoRing) -> (GrpAbFinGen, GrpAbFinGenToAbsOrdQuoRingMultMap)

Returns the group $Q^\times$ and a map from this group to $Q$.
"""
function _multgrp(Q::NfOrdQuoRing, save_tame_wild::Bool = false; method = nothing)

  fac = factor(Q)

  if isempty(fac)
    disc_log = function(x::NfOrdQuoRingElem)
      return fmpz[]
    end
    GtoQ = GrpAbFinGenToAbsOrdQuoRingMultMap(Q, elem_type(Q)[], fmpz[], disc_log)
    return domain(GtoQ), GtoQ
  end

  prime_powers = Vector{NfOrdIdl}()
  groups = Vector{GrpAbFinGen}()
  maps = Vector{GrpAbFinGenToNfOrdQuoRingMultMap}()
  tame_ind = Tuple{NfOrdIdl, Int}[]
  ind = 1
  for (p, vp) in fac
    pvp = p^vp
    G1, mG1 = _multgrp_mod_pv(p, vp, pvp; method = method)
    push!(tame_ind, (p, ind))
    ind += ngens(G1)
    push!(prime_powers, pvp)
    push!(groups, G1)
    push!(maps, mG1)
  end

  G, GtoQ = _direct_product(groups, maps, prime_powers, Q, save_tame_wild)
  S, StoG, StoQ = snf(G, GtoQ)

  if save_tame_wild
    for s = 1:length(tame_ind)
      StoQ.tame[tame_ind[s][1]].disc_log = StoG\(G[tame_ind[s][2]]) #Relies on the ordering tame\wild in the construction!
    end
  end

  return S, StoQ
end

################################################################################
#
#  Compute Multiplicative Group For Prime Powers
#
################################################################################

@doc Markdown.doc"""
    _multgrp_mod_pv(p::NfOrdIdl, v::Int, pv::NfOrdIdl) -> (GrpAbFinGen, GrpAbFinGenToAbsOrdQuoRingMultMap)

Given a prime ideal $p$ in a maximal order $\mathcal O$, an integer $v > 0$ and
$pv = p^v$, the function returns the group $(\mathcal O/p^v)^\times$ and a map
from this group to $O/p^v$.
"""
function _multgrp_mod_pv(p::NfOrdIdl, v::Int, pv::NfOrdIdl; method=nothing)
  @hassert :NfOrdQuoRing 2 is_prime(p)
  @assert v >= 1
  pnumv = minimum(p, copy = false)^v # to speed up the exponentiation in the GrpAbFinGenToNfAbsOrdMaps
  G1, G1toO = _multgrp_mod_p(p, pnumv)
  Q, OtoQ = quo(order(p), pv)
  tame_part = Dict{NfAbsOrdIdl, GrpAbFinGenToNfAbsOrdMap}()
  wild_part = Dict{NfAbsOrdIdl, GrpAbFinGenToNfAbsOrdMap}()
  if v == 1
    G1toO.disc_log = G1[1]
    tame_part[p] = G1toO
    function disc_log(x::NfOrdQuoRingElem)
      return G1toO.discrete_logarithm((OtoQ\x))
    end
    GtoQ = GrpAbFinGenToAbsOrdQuoRingMultMap(G1, Q, map(OtoQ, G1toO.generators), disc_log)
  else
    G2, G2toO = _1_plus_p_mod_1_plus_pv(p, v, pv, pnumv; method = method)
    @assert is_snf(G2)
    wild_part[p] = G2toO
    gen1 = G1toO(G1[1])

    rel1 = exponent(G1)
    # G2.snf[end] is the order of the biggest cyclic subgroup of G2
    gen1_obcs = powermod(gen1, G2.snf[end], pnumv)
    gens = map(OtoQ, append!([gen1_obcs], G2toO.generators))

    G1toO.generators[1] = gen1_obcs

    G = direct_product(G1, G2, task = :none)

    G1toO.disc_log = G[1]
    tame_part[p] = G1toO

    obcs_inv = gcdx(G2.snf[end], rel1)[2]
    function disc_log2(x::NfOrdQuoRingElem)
      y = OtoQ\x
      r = mod((G1toO.discrete_logarithm(y))[1]*obcs_inv, rel1)
      y *= powermod(gen1_obcs, mod(-r, rel1), pnumv)
      return append!([r], G2toO.discrete_logarithm(y))
    end

    GtoQ = GrpAbFinGenToAbsOrdQuoRingMultMap(G, Q, gens, disc_log2)
  end
  GtoQ.tame = tame_part
  GtoQ.wild = wild_part
  return domain(GtoQ), GtoQ
end

################################################################################
#
#  Compute Multiplicative Group For Primes
#
################################################################################

# Compute (O_K/p)*
function _multgrp_mod_p(p::NfOrdIdl, pnumv::fmpz = fmpz(0))
  @hassert :NfOrdQuoRing 2 is_prime(p)
  O = order(p)
  n = norm(p) - 1
  gen = _primitive_element_mod_p(p)
  factor_n = factor(n)
  Q, mQ = ResidueField(O, p)
  gen_quo = mQ(gen)
  big_step_cache = Dict{fmpz, Dict{typeof(gen_quo), fmpz}}()
  local discrete_logarithm
  let mQ = mQ, n = n, Q = Q, gen_quo = gen_quo, factor_n = factor_n, big_step_cache = big_step_cache
    function discrete_logarithm(x::NfOrdElem)
      y = mQ(x)
      if isone(y)
        return fmpz[0]
      elseif y == Q(-1) && iszero(mod(n, 2))
        return fmpz[divexact(n, 2)]
      end
      if n < 11
        res = 1
        el = gen_quo
        while el != y
          el *= gen_quo
          res += 1
        end
        return fmpz[res]
      else
        return fmpz[pohlig_hellman(gen_quo,n,y;factor_n=factor_n, big_step_cache = big_step_cache)]
      end
    end
  end
  map = GrpAbFinGenToNfAbsOrdMap(O, [ gen ], [ n ], discrete_logarithm, pnumv)
  return domain(map), map
end

function _primitive_element_mod_p(p::NfOrdIdl)
  @hassert :NfOrdQuoRing 2 is_prime(p)
  O = order(p)
  Q, Q_map = quo(O,p)
  n = norm(p) - 1
  primefactors_n = collect(keys(factor(n).fac))
  while true
    x = rand(Q)
    x == 0 && continue
    order_too_small = false
    for l in primefactors_n
      if isone(x^div(n, l))
        order_too_small = true
        break
      end
    end
    order_too_small || return Q_map\x
  end
end

################################################################################
#
# Computation of (1+p)/(1+p^v)
#
################################################################################

# Compute (1+p)/(1+p^v)
function _1_plus_p_mod_1_plus_pv(p::NfOrdIdl, v::Int, pv::NfOrdIdl, pnumv::fmpz = fmpz(0); method=nothing)
  @hassert :NfOrdQuoRing 2 is_prime(p)
  @assert v >= 1

  if method == :quadratic
    gens, rels, disc_log = _iterative_method(p,v;base_method=:quadratic,use_p_adic=false)
  elseif method == :artin_hasse
    gens, rels, disc_log = _iterative_method(p,v;base_method=:artin_hasse,use_p_adic=false)
  elseif method == :p_adic
    gens, rels, disc_log = _iterative_method(p,v;use_p_adic=true)
  else
    gens, rels, disc_log = _iterative_method(p,v)
  end

  @assert size(rels) == (length(gens),length(gens))
  toO = GrpAbFinGenToNfAbsOrdMap(order(p), gens, rels, disc_log, pnumv)::GrpAbFinGenToAbsOrdMap{NfOrd, NfOrdElem}
  Q, OtoQ = quo(order(p), pv)
  S, mS, StoO = snf(toO, OtoQ)
  return S, StoO::GrpAbFinGenToAbsOrdMap{NfOrd, NfOrdElem}
end

################################################################################
#
#  Iterative Method for (1+p^u)/(1+p^v)
#
################################################################################

function _iterative_method(p::NfOrdIdl, v::Int; base_method=nothing, use_p_adic=true)
  return _iterative_method(p, 1, v ; base_method = base_method, use_p_adic = use_p_adic)
end

function _iterative_method(p::NfOrdIdl, u::Int, v::Int; base_method = nothing, use_p_adic = true)
  @assert v >= u >= 1
  pnum = minimum(p)
  if use_p_adic
    e = ramification_index(p)
    k0 = 1 + Int(div(fmpz(e),(pnum-1)))::Int
  end
  g = Vector{NfOrdElem}()
  M = zero_matrix(FlintZZ, 0, 0)
  dlogs = Vector{Function}()

  l = u
  pl = (p^l)::NfOrdIdl

  while l != v
    k = l
    pk = pl
    if use_p_adic && k>=k0
      l = v
      pl = p^l
      h, N, disc_log = _p_adic_method(p, k, l; pu = pk, pv = pl)::Tuple{Vector{NfOrdElem}, Vector{fmpz}, Function}
    elseif base_method == :quadratic
      l = min(2*k, v)
      pl = p^l
      h, N, disc_log = _quadratic_method(p, k, l; pu = pk, pv = pl)::Tuple{Vector{NfOrdElem}, Vector{fmpz}, Function}
    else
      l = Int(min(pnum*k, v))
      pl = p^l
      h, N, disc_log = _artin_hasse_method(p, k, l; pu = pk, pv = pl)::Tuple{Vector{NfOrdElem}, Vector{fmpz}, Function}
    end
    g, M = _expand(g, M, h, N, disc_log, pl)
    push!(dlogs, disc_log)
  end

  Q = NfOrdQuoRing(order(pl), pl)
  local discrete_logarithm
  let Q = Q, dlogs = dlogs, pl = pl
    function discrete_logarithm(b::NfOrdElem)
      b1 = Q(b)
      a = fmpz[]
      k1 = 1
      for i in 1:length(dlogs)
        a_ = dlogs[i](b1.elem)
        prod = one(Q)
        for j in 1:length(a_)
          if !iszero(a_[j])
            mul!(prod, prod, Q(g[k1])^a_[j])
          end
          k1 += 1
        end
        append!(a, a_)
        b1 = divexact(b1, prod)
      end
      return a
    end
  end

  return g, M, discrete_logarithm
end

#
# Given generators and relations for groups of two consecutives steps, this function computes
# generators and relations for the product
#
function _expand(g::Vector{NfOrdElem}, M::fmpz_mat, h::Vector{NfOrdElem}, N::Vector{fmpz}, disc_log::Function, pl::NfOrdIdl)
  if isempty(g)
    M1 = zero_matrix(FlintZZ, length(N), length(N))
    for i = 1:length(N)
      M1[i, i] = N[i]
    end
    return h, M1
  end
  isempty(h) && return g,M
  O = order(pl)
  Q , mQ = quo(O, pl)
  Z = zero_matrix(FlintZZ,nrows(M)+length(N),ncols(M)+length(N))
  for i = 1:nrows(M)
    for j = i:ncols(M)
      Z[i, j] = M[i, j]
    end
  end
  for i=1:length(N)
    Z[i+nrows(M),i+nrows(M)] = N[i]
  end
  for i in 1:nrows(M)
    el = one(Q)
    for j = 1:ncols(M)
      if !iszero(M[i, j])
        mul!(el, el, Q(g[j])^M[i, j])
      end
    end
    alpha = disc_log(el.elem)
    for j in 1:length(N)
      Z[i,j+ncols(M)] = -alpha[j]
    end
  end
  append!(g,h)

  return g, Z
end

#
#  This function returns a set of generators with the corresponding relations and disclog
#

function _pu_mod_pv(pu::NfOrdIdl, pv::NfOrdIdl)

  O=order(pu)
  b=basis(pu)
  N = basis_matrix(pv, copy = false)*basis_mat_inv(pu, copy = false)
  @assert isone(N.den)
  G = abelian_group(N.num)
  S, mS=snf(G)

  #Generators
  gens=Vector{NfOrdElem}(undef, ngens(S))
  for i=1:ngens(S)
    x=mS(S[i])
    gens[i]= zero(O)
    for j=1:ngens(G)
      add!(gens[i], gens[i], mod(x[j], S.snf[end])*b[j])
    end
  end

  #Disclog
  M = basis_mat_inv(pu, copy = false)*mS.imap
  x_fakemat2 = FakeFmpqMat(zero_matrix(FlintZZ, 1, ncols(M)), fmpz(1))
  local disclog
  let M = M, O = O, S = S, x_fakemat2 = x_fakemat2
    function disclog(x::NfOrdElem)
      x_fakemat = FakeFmpqMat(matrix(FlintZZ, 1, degree(O), coordinates(x, copy = false)))
      mul!(x_fakemat2, x_fakemat, M)
      #@assert x_fakemat2 == x_fakemat * M
      denominator(x_fakemat2) != 1 && error("Element is in the ideal")
      res = Vector{fmpz}(undef, ncols(x_fakemat2))
      for i = 1:length(res)
        res[i] = x_fakemat2.num[1, i]#mod(x_fakemat.num[1, i], S.snf[end])
      end
      return res
    end
  end
  return gens, S.snf, disclog

end

# Let p be a prime ideal above a prime number pnum. Let e = v_p(pnum) be
# its ramification index. If b > a >= e/(pnum-1) this function computes
# the structure of (1+p^a)/(1+p^b) as an abelian group.
function _1_plus_pa_mod_1_plus_pb_structure(p::NfOrdIdl,a,b)
  b > a >= 1 || return false, nothing
  @hassert :NfOrdQuoRing 2 is_prime(p)
  O = order(p)
  pnum = minimum(p)
  e = ramification_index(p)
  k0 = 1 + div(fmpz(e),(pnum-1))
  a >= k0 || return false, nothing
  Q = NfOrdQuoRing(O,p^(b-a))
  return true, group_structure(Q)
end

################################################################################
#
# Quadratic Method for (1+p^u)/(1+p^v)
#
################################################################################

# Compute generators, a relation matrix and a function to compute discrete
# logarithms for (1+p^u)/(1+p^v), where 2*u >= v >= u >= 1
function _quadratic_method(p::NfOrdIdl, u::Int, v::Int; pu::NfOrdIdl=p^u, pv::NfOrdIdl=p^v)
  @hassert :NfOrdQuoRing 2 is_prime(p)
  @assert 2*u >= v >= u >= 1
  g, M, dlog = _pu_mod_pv(pu,pv)
  map!(x -> x + 1, g, g)
  function discrete_logarithm(x::NfOrdElem)
    res = dlog(mod(x-1,pv))
    for i = 1:length(res)
      res[i] = mod(res[i], M[end])
    end
    return res
  end
  return g, M, discrete_logarithm
end

################################################################################
#
# Artin-Hasse Method for (1+p^u)/(1+p^v)
#
################################################################################

# Compute generators, a relation matrix and a function to compute discrete
# logarithms for (1+p^u)/(1+p^v), where p is a prime ideal over pnum
# and pnum*u >= v >= u >= 1
function _artin_hasse_method(p::NfOrdIdl, u::Int, v::Int; pu::NfOrdIdl=p^u, pv::NfOrdIdl=p^v)
  @hassert :NfOrdQuoRing 2 is_prime(p)
  pnum = minimum(p)
  @assert pnum*u >= v >= u >= 1
  @assert fmpz(v) <= pnum*fmpz(u)
  @assert has_minimum(pv)
  Q, mQ = quo(order(p), pv)
  g, M, dlog = _pu_mod_pv(pu,pv)
  map!(x -> artin_hasse_exp(Q(x), pnum), g, g)
  local discrete_logarithm
  let Q = Q, pnum = pnum, dlog = dlog
    function discrete_logarithm(x::NfOrdElem)
      res = dlog(artin_hasse_log(Q(x), pnum))
      for i = 1:length(res)
        res[i] = mod(res[i], M[end])
      end
      return res
    end
  end
  return g, M, discrete_logarithm
end

function artin_hasse_exp(x::NfOrdQuoRingElem, pnum::fmpz)
  Q = parent(x)
  s = Q(1)
  fac_i = fmpz(1)
  t = Q(1)
  m = minimum(Q.ideal)
  for i=1:pnum-1
    mul!(t, t, x)
    if iszero(t)
      break
    end
    mul!(fac_i, fac_i, i)
    mod!(fac_i, fac_i, m)
    invfac_i = invmod(fac_i, m)
    add!(s, s, invfac_i*t)
  end
  return s.elem
end

function artin_hasse_log(y::NfOrdQuoRingElem, pnum::fmpz)
  Q = parent(y)
  x = y - Q(1)
  if iszero(x)
    return x.elem
  end
  m = minimum(Q.ideal)
  s = Q(0)
  t = Q(1)
  for i in 1:pnum-1
    mul!(t, t, x)
    if iszero(t)
      break
    end
    invi = invmod(fmpz(i), m)
    if iseven(i)
      sub!(s, s, invi*t)
    else
      add!(s, s, invi*t)
    end
  end
  return s.elem
end

################################################################################
#
# p-Adic Method for (1+p^u)/(1+p^v)
#
################################################################################

# Compute generators, a relation matrix and a function to compute discrete
# logarithms for (1+p)/(1+p^v) if 1 >= k0, where p is a prime ideal over pnum,
# e the p-adic valuation of pnum, and k0 = 1 + div(e,pnum-1)
function _p_adic_method(p::NfOrdIdl, v::Int; pv=p^v)
  return _p_adic_method(p,1,v)
end

# Compute generators, a relation matrix and a function to compute discrete
# logarithms for (1+p^u)/(1+p^v) if u >= k0, where p is a prime ideal over pnum,
# e the p-adic valuation of pnum, and k0 = 1 + div(e,pnum-1)
function _p_adic_method(p::NfOrdIdl, u::Int, v::Int; pu::NfOrdIdl=p^u, pv::NfOrdIdl=p^v)
  @assert v > u >= 1
  e = ramification_index(p) #ramification index
  k0 = 1 + div(fmpz(e),(minimum(p)-1))
  @assert u >= k0
  g, M, dlog = _pu_mod_pv(pu, pv)
  O = order(p)
  Q = NfOrdQuoRing(O, pv)
  map!(x -> p_adic_exp(Q, p, v, x), g, g)
  powers = Dict{Int, nf_elem}()
  local discrete_logarithm
  let Q = Q, p = p, v = v, dlog = dlog, powers = powers
    function discrete_logarithm(b::NfOrdElem)
      res = dlog(p_adic_log(Q, p, v, b, powers))
      for i = 1:length(res)
        res[i] = mod(res[i], M[end])
      end
      return res
    end
  end
  return g, M, discrete_logarithm
end

function _divexact(x::AbsOrdQuoRingElem, y::fmpz)
  Q = parent(x)
  I = ideal(Q)
  OK = order(I)
  m = minimum(I, copy = false)
  s = invmod(y, m)
  res = s*x
  #@assert res == divexact(x, Q(y))
  return res
end

function p_adic_exp(Q::NfOrdQuoRing, p::NfOrdIdl, v::Int, x::NfOrdElem)
  O = parent(x)
  iszero(x) && return one(O)
  pnum = p.minimum
  e = ramification_index(p)
  val_p_x = valuation(x, p)
  max_i = floor(Int, v / (val_p_x - (e/(Float64(pnum)-1)))) + 1
  val_p_maximum = Int(max_i*val_p_x) + 1
  Q_ = NfOrdQuoRing(O, p^val_p_maximum)
  x1 = Q_(x)
  s = one(Q)
  inc = one(Q_)
  val_p_xi = 0
  val_p_fac_i = 0
  i_old = 0
  for i in 1:max_i
    val_pnum_i = valuation(fmpz(i), pnum)
    val_p_i = val_pnum_i * e
    val_p_fac_i += val_p_i
    val_p_xi += val_p_x
    val_p_xi - val_p_fac_i >= v && continue
    @hassert :NfOrdQuoRing 1 val_p_xi - val_p_fac_i>=0
    @hassert :NfOrdQuoRing 1 val_p_xi< val_p_maximum
    @hassert :NfOrdQuoRing 1 val_p_fac_i< val_p_maximum
    i_prod = prod((i_old+1):i)
    deltax = inc*x1^(i-i_old)
    @hassert :NfOrdQuoRing 1 !iszero(deltax)
    if isone(gcd(i_prod, pnum))
      inc = _divexact(deltax, fmpz(i_prod))
    else
      inc = divexact(deltax, Q_(i_prod))
    end
    add!(s, s, Q(inc.elem))
    i_old = i
  end
  return s.elem
end

function p_adic_log(Q::NfOrdQuoRing, p::NfOrdIdl, v::Int, y::NfOrdElem, powers::Dict{Int, nf_elem} = Dict{Int, nf_elem}())
  O = parent(y)
  isone(y) && return zero(O)
  pnum = minimum(p)
  x = y - 1
  val_p_x = valuation(x, p)
  s = zero(Q)
  xi = one(O)
  i_old = 0
  val_p_xi = 0
  anti_uni = anti_uniformizer(p)
  e = ramification_index(p)
  #we have to find a bound for this.
  # it is enough to find the minimum l such that
  # pnum^l v_p(x)>= v+el
  l = 1
  left = pnum*val_p_x
  right=v+e
  while left < right
    left*=pnum
    right+=e
    l+=1
  end
  bound2=pnum^l-pnum
  bound1 = div(v, val_p_x)
  leftlim = bound1 + pnum - mod(fmpz(bound1), pnum)
  for i in 1:bound1
    val_pnum_i = valuation(i, pnum)
    val_p_i = val_pnum_i * e
    val_p_xi += val_p_x
    val_p_xi - val_p_i >= v && continue
    mul!(xi, xi, x^(i-i_old))
    if iszero(val_pnum_i)
      inc = _divexact(Q(xi), fmpz(i))
    else
      if haskey(powers, val_p_i)
        el = powers[val_p_i]
      else
        el = _powermod(anti_uni, val_p_i, pnum^(val_p_i+1))
        powers[val_p_i] = el
      end
      numer = O(xi.elem_in_nf*el, false)
      denom = O(i*el, false)
      inc = divexact(Q(numer),Q(denom))
    end
    if isodd(i)
      add!(s, s, inc)
    else
      sub!(s, s, inc)
    end
    i_old = i
  end
  for i in leftlim:pnum:bound2
    val_pnum_i = valuation(i, pnum)
    val_p_i = val_pnum_i * e
    val_p_xi += val_p_x
    val_p_xi - val_p_i >= v && continue
    mul!(xi, xi, x^(i-i_old))
    if iszero(val_pnum_i)
      inc = _divexact(Q(xi), fmpz(i))
    else
      if haskey(powers, val_p_i)
        el = powers[val_p_i]
      else
        el = _powermod(anti_uni, val_p_i, pnum^(val_p_i+1))
        powers[val_p_i] = el
      end
      numer = O(xi.elem_in_nf*el, false)
      denom = O(i*el, false)
      inc = divexact(Q(numer),Q(denom))
    end
    if isodd(i)
      add!(s, s, inc)
    else
      sub!(s, s, inc)
    end
    i_old = i
  end
  return s.elem
end

################################################################################
#
#  Discrete Logarithm In Cyclic Groups
#
################################################################################
# TODO compare with implementations in UnitsModM.jl

@doc Markdown.doc"""
    baby_step_giant_step(g, n, h) -> fmpz
    baby_step_giant_step(g, n, h, cache::Dict) -> fmpz

Computes the discrete logarithm $x$ such that $h = g^x$.

$g$ is a generator of order less than or equal to $n$
and $h$ has to be generated by $g$.
If a dictionary `cache` is supplied, it will be used to store the result
of the first step. This allows to speed up subsequent calls with
the same $g$ and $n$.
"""
function baby_step_giant_step(g, n, h, cache::Dict)
  @assert typeof(g) == typeof(h)
  m = fmpz(root(n, 2)+1)
  if isempty(cache)
    it = g^0
    for j in 0:m
      cache[it] = j
      it *= g
    end
  end
  b = g^(-m)
  y = h
  for i in 0:m-1
    if haskey(cache, y)
      return fmpz(mod(i*m + cache[y], n))
    else
      y *= b
    end
  end
  error("Couldn't find discrete logarithm")
end

function baby_step_giant_step(gen, n, a)
  cache = Dict{typeof(gen), BigInt}()
  return baby_step_giant_step(gen, n, a, cache)
end

@doc Markdown.doc"""
    pohlig_hellman(g, n, h; factor_n=factor(n)) -> fmpz

Computes the discrete logarithm $x$ such that $h = g^x$.

$g$ is a generator of order $n$ and $h$ has to be generated by $g$.
The factorisation of $n$ can be supplied via `factor_n` if
it is already known.
"""
function pohlig_hellman(g, n, h; factor_n=factor(n), big_step_cache = Dict())
  @assert typeof(g) == typeof(h)
  n == 1 && return fmpz(0)
  results = Vector{fmpz}()
  prime_powers = Vector{fmpz}()
  for (p,v) in factor_n
    pv = p^v
    r = div(n,pv)
    if !haskey(big_step_cache, p)
      big_step_cache[p] = Dict{typeof(g), fmpz}()
    end
    c = _pohlig_hellman_prime_power(g^r,p,v,h^r, big_step_cache = big_step_cache[p])
    push!(results, fmpz(c))
    push!(prime_powers, fmpz(pv))
  end
  if length(results) == 1
    return results[1]
  end
  return crt(results, prime_powers)
end

function _pohlig_hellman_prime_power(g, p, v, h; big_step_cache::Dict)
  p_i = 1
  p_v_min_i_min_1 = p^(v-1)
  g_ = g^(p^(v-1))
  a = baby_step_giant_step(g_,p,h^(p^(v-1)), big_step_cache)
  h *= g^-a
  for i in 1:v-1
    p_i *= p
    p_v_min_i_min_1 = div(p_v_min_i_min_1,p)
    ai = baby_step_giant_step(g_, p , h^p_v_min_i_min_1, big_step_cache)
    ai_p_i = ai * p_i
    a += ai_p_i
    h *= g^(-ai_p_i)
  end
  return a
end

#################################################################################
#
#  p-part of the multiplicative group
#
#################################################################################

function _prime_part_multgrp_mod_p(p::NfOrdIdl, prime::Int)
  @hassert :NfOrdQuoRing 2 is_prime(p)
  O = order(p)
  Q, mQ = ResidueField(O,p)

  n = norm(p) - 1
  s=valuation(n,prime)
  powerp=prime^s
  m=divexact(n,powerp)

  powm=divexact(powerp,prime)
  found=false
  g=Q(1)
  while found==false
    g = rand(Q)
    if g != Q(0)
      g=g^m
      if g^powm != Q(1)
        found=true
      end
    end
  end
  inv=gcdx(m,fmpz(powerp))[2]

  function disclog(x::NfOrdElem)
    t=mQ(x)^m
    if powerp<10
      w=1
      el=g
      while el!=t
        w+=1
        el*=g
      end
      return fmpz[w*inv]
    else
      return fmpz[_pohlig_hellman_prime_power(g,prime,s,t)*inv]
    end
  end

  map = GrpAbFinGenToNfAbsOrdMap(O, [ mQ\g ], [ fmpz(powerp) ], disclog)
  return domain(map), map
end

function _mult_grp(Q::NfOrdQuoRing, p::Integer)
  O = Q.base_ring
  OtoQ = NfOrdQuoMap(O, Q)

  tame_part = Dict{NfOrdIdl, GrpAbFinGenToNfAbsOrdMap}()
  wild_part = Dict{NfOrdIdl, GrpAbFinGenToNfAbsOrdMap}()

  fac = factor(Q)
  y1 = Dict{NfOrdIdl, Int}()
  y2 = Dict{NfOrdIdl, Int}()
  for (q, e) in fac
    if divisible(norm(q) - 1, p)
      y1[q] = Int(1)
    else
      if divisible(norm(q), p) && e >= 2
        y2[q] = Int(e)
      end
    end
  end

  if isempty(y1) && isempty(y2)
    G = abelian_group(fmpz[])
    disc_log = function(x::NfOrdQuoRingElem)
      return fmpz[]
    end
    return G, GrpAbFinGenToAbsOrdQuoRingMultMap(G, Q, elem_type(Q)[], disc_log), y1
  end

  ideals_and_maps = Dict{NfOrdIdl, Vector{GrpAbFinGenToNfAbsOrdMap}}()
  for (q, e) in fac
    qe = q^e
    maps = GrpAbFinGenToNfAbsOrdMap[]

    if haskey(y1, q)
      G1, G1toO = _prime_part_multgrp_mod_p(q, p)
      push!(maps, G1toO)
      tame_part[q] = G1toO
    end

    if haskey(y2, q)
      G2, G2toO = _1_plus_p_mod_1_plus_pv(q, y2[q], q^y2[q])

      nq = norm(q) - 1

      @assert is_snf(G2)
      obcs = G2.snf[end] # order of the biggest cyclic subgroup
      obcs_inv = gcdx(nq, obcs)[2]

      disc_log2 = function(x::NfOrdElem)
        y = Q(x)^Int(nq)
        z = G2toO.discrete_logarithm(y.elem)
        for i = 1:length(z)
          z[i] *= obcs_inv
        end
        return z
      end

      G2toO2 = GrpAbFinGenToNfAbsOrdMap(G2, O, G2toO.generators, disc_log2)
      push!(maps, G2toO2)
      wild_part[q] = G2toO2
    end
    ideals_and_maps[qe] = maps
  end
  G, GtoQ = _direct_product!(ideals_and_maps, Q) # This also changes tame_part and wild_part!
  GtoQ.tame = tame_part
  GtoQ.wild = wild_part

  S, StoG, StoQ = snf(G, GtoQ)

  return S, StoQ, merge(y1, y2)
end

####################################################################################
#
#  multiplicative group mod n
#
####################################################################################

function _find_gen(Q::FqFiniteField, powm::Vector{fmpz}, m::fmpz)
  found = false
  g = one(Q)
  while !found
    g = rand(Q)
    !iszero(g) || continue
    g = g^m
    found = true
    for i=1:length(powm)
      if isone(g^powm[i])
        found = false
        break
      end
    end
  end
  return g
end

function _n_part_multgrp_mod_p(p::NfOrdIdl, n::Int)
  @hassert :NfOrdQuoRing 2 is_prime(p)
  O = order(p)
  Q, mQ = ResidueField(O, p)

  np = norm(p) - fmpz(1)
  @assert !isone(gcd(fmpz(n), np))
  npart, m = ppio(np, fmpz(n))
  k = gcd(npart, fmpz(n))
  fac = factor(k)
  powm = fmpz[divexact(npart, x) for x in keys(fac.fac)]

  #
  #  We search for a random element with the right order
  #
  g = _find_gen(Q, powm, m)
  inv = gcdx(m, npart)[2]
  quot = divexact(npart, k)

  w = g^quot
  local disclog
  let m = m, quot = quot, k = k, w = w, inv = inv
    function disclog(x::NfOrdElem)
      t = mQ(x)^(m*quot)
      if iszero(t)
        error("Not coprime!")
      end
      if isone(t)
        return fmpz[fmpz(0)]
      end
      if k < 20
        s = 1
        el = deepcopy(w)
        while el != t
          s += 1
          mul!(el, el, w)
        end
        return fmpz[mod(fmpz(s)*inv, k)]
      else
        return fmpz[pohlig_hellman(w, k, t)*inv]
      end
    end
  end
  G = abelian_group([k])
  gens = Vector{NfOrdElem}(undef, 1)
  gens[1] = preimage(mQ, g)
  map = GrpAbFinGenToNfAbsOrdMap(G, O, gens, disclog)::GrpAbFinGenToAbsOrdMap{NfOrd, NfOrdElem}
  return G, map
end

#Cohen, Advanced topics in computational number theory, 4.5 exercise 20
function bound_exp_mult_grp(P::NfOrdIdl, n::Int)
  @assert is_prime(P)
  e = ramification_index(P)
  f = degree(P)
  p = minimum(P, copy = false)
  s = valuation(n, p)
  return Int((p^s)-1+s*p^s)
end

function _mult_grp_mod_n(Q::NfOrdQuoRing, y1::Dict{NfOrdIdl, Int}, y2::Dict{NfOrdIdl, Int}, n::Integer)

  O = Q.base_ring
  idQ = Q.ideal
  OtoQ = NfOrdQuoMap(O, Q)
  fac = factor(Q)

  tame_part = Dict{NfOrdIdl, GrpAbFinGenToNfAbsOrdMap}()
  wild_part = Dict{NfOrdIdl, GrpAbFinGenToNfAbsOrdMap}()

  if isempty(y1) && isempty(y2)
    G = abelian_group(fmpz[])
    disc_log = function(x::NfOrdQuoRingElem)
      return fmpz[]
    end
    return G, GrpAbFinGenToAbsOrdQuoRingMultMap(G, Q, elem_type(Q)[], disc_log), tame_part, wild_part
  end

  ideals_and_maps = Tuple{NfOrdIdl, Vector{GrpAbFinGenToNfAbsOrdMap}}[]
  tame_ind = Tuple{NfOrdIdl, Int}[]
  i = 0
  for (q, e) in fac
    qe = q^e
    maps = GrpAbFinGenToNfAbsOrdMap[]

    if haskey(y1, q)
      G1, G1toO = _n_part_multgrp_mod_p(q, n)
      push!(maps, G1toO)
      tame_part[q] = G1toO
      i += 1
      push!(tame_ind, (q, i))
    end

    if haskey(y2, q)
      @assert y2[q] >= 2
      #exp_q = min(y2[q], bound_exp_mult_grp(q, n))
      G2, G2toO = _1_plus_p_mod_1_plus_pv(q, y2[q], qe)
      if haskey(y1, q)
        e = exponent(G2)
        uncom = Int(ppio(fmpz(n), e)[2])
        e1 = e
        while !isone(mod(e1, uncom))
          e1 *= e
        end
        tame_part[q].generators[1] = powermod(tame_part[q].generators[1], e1, minimum(idQ, copy = false))
      end

      i += ngens(G2)
      nq = norm(q) - 1

      @assert is_snf(G2)
      obcs = G2.snf[end] # order of the biggest cyclic subgroup
      obcs_inv = gcdx(nq, obcs)[2]

      local disc_log2
      let Q = Q, nq = nq, G2toO = G2toO
        function disc_log2(x::NfOrdElem)
          y = Q(x)^nq
          z = G2toO.discrete_logarithm(y.elem)
          for i = 1:length(z)
            z[i] *= obcs_inv
          end
          return z
        end
      end

      G2toO2 = GrpAbFinGenToNfAbsOrdMap(G2, O, G2toO.generators, disc_log2)::GrpAbFinGenToAbsOrdMap{NfOrd, NfOrdElem}
      push!(maps, G2toO2)
      wild_part[q] = G2toO2
    end
    push!(ideals_and_maps, (qe, maps))
  end
  G, GtoQ = _direct_product!(ideals_and_maps, Q) # This also changes tame_part and wild_part!

  S, StoG, StoQ = snf(G, GtoQ)

  for s = 1:length(tame_ind)
    tame_part[tame_ind[s][1]].disc_log = StoG\(G[tame_ind[s][2]])
  end
  StoQ.tame = tame_part
  StoQ.wild = wild_part
  return S, StoQ
end

################################################################################
#
#  Functions for the non-maximal case
#
################################################################################

# Computes generators and relations for (O_P/AO_P)^\times, where Pm = P^m, Q = O/A.
function _multgrp_Op_aOp(Q::NfOrdQuoRing, P::NfOrdIdl, m::Int, Pm::NfOrdIdl)
  A = ideal(Q)
  O = order(A)

  # Compute generators for (1 + (A + P^m))/(1 + P^m)
  I = A + Pm
  I2 = I^2
  S = Set{NfOrdElem}()
  j = 0
  k = floor(Int, log(2, m))
  while j <= k
    if j > 0
      I = I2
      I2 = I2^2
    end
    for b in basis(I)
      push!(S, mod(b, I2))
    end
    j += 1
  end
  O1 = order(A)(1)
  gens = [ g + O1 for g in S ]

  # Compute (O/P^m)^\times
  G, GtoQQ = _multgrp_mod_pv(P, m, Pm, method = :quadratic)
  QQ = codomain(GtoQQ)
  OtoQQ = NfOrdQuoMap(O, QQ)

  H, GtoH = quo(G, [ GtoQQ\OtoQQ(g) for g in gens])
  @assert GtoH.map == identity_matrix(FlintZZ, ngens(G))
  HtoQQ = GrpAbFinGenToAbsOrdQuoRingMultMap(H, QQ, GtoQQ.generators, GtoQQ.discrete_logarithm)

  return H, HtoQQ
end

function _multgrp_non_maximal(Q::NfOrdQuoRing)
  A = ideal(Q)
  O = order(A)

  # Determine the prime ideals containing A
  OO = maximal_order(nf(O))
  aOO = extend(A, OO)
  fac = factor(aOO)
  S = Set{NfOrdIdl}()
  for (p, e) in fac
    q = contract(p, O)
    q.is_prime = p.is_prime
    push!(S, q)
  end
  prime_ideals = [ p for p in S ]

  # Compute (upper bounds for) exponents m such that A*O_P \supseteq P^m*O_P
  m = Vector{Int}(undef, length(prime_ideals))
  for i = 1:length(m)
    P = prime_ideals[i]
    p = minimum(P)
    x = valuation(det(basis_matrix(A, copy = false)), p)
    y = valuation(det(basis_matrix(P, copy = false)), p)
    m[i] = Int(ceil(x/y))
  end

  # Compute the groups (O_P/AO_P)^\times
  Pm = [ prime_ideals[i]^m[i] for i in 1:length(prime_ideals)]
  groups = Vector{GrpAbFinGen}(undef, length(prime_ideals))
  maps = Vector{GrpAbFinGenToNfOrdQuoRingMultMap}(undef, length(prime_ideals))
  for i= 1:length(prime_ideals)
    H, HtoQQ = _multgrp_Op_aOp(Q, prime_ideals[i], m[i], Pm[i])
    groups[i] = H
    maps[i] = HtoQQ
  end

  G, GtoQ = _direct_product(groups, maps, Pm, Q)

  S, StoG, StoQ = snf(G, GtoQ)
  return S, StoQ
end

################################################################################
#
#  SNF and direct product (which transform the map too)
#
################################################################################

# For an element x of elements[i] this computes an element y with
# x \equiv y mod ideals[i] and x \equiv 1 mod ideals[j] for all j not equal i.
function make_coprime(elements::Vector{Vector{S}}, ideals::Vector{T}) where { S <: Union{ NfAbsOrdElem, AlgAssAbsOrdElem }, T <: Union{ NfAbsOrdIdl, AlgAssAbsOrdIdl } }
  @assert !isempty(ideals)
  @assert length(elements) == length(ideals)

  n = length(ideals)
  if n == 1
    return elements[1]
  end

  products = _compute_products_for_make_coprime(ideals)

  One = one(order(ideals[1]))
  result = Vector{S}()
  for i = 1:n
    for j = 1:length(elements[i])
      t = crt(elements[i][j], ideals[i], One, products[i])
      push!(result, t)
    end
  end
  return result
end

# Build the products \prod_{j\neq i} ideals[j] for all i
function _compute_products_for_make_coprime(ideals::Vector{T}) where { T <: Union{ NfAbsOrdIdl, AlgAssAbsOrdIdl } }
  n = length(ideals)
  @assert n >= 2

  left_halves = Vector{T}(undef, n - 1)
  right_halves = Vector{T}(undef, n - 1)
  left_halves[1] = ideals[1]
  right_halves[1] = ideals[n]
  for i = 2:(n - 1)
    left_halves[i] = left_halves[i - 1]*ideals[i]
    right_halves[i] = right_halves[i - 1]*ideals[n - i + 1]
  end
  products = Vector{T}(undef, n)
  products[1] = right_halves[n - 1]
  products[n] = left_halves[n - 1]
  for i = 2:(n - 1)
    products[i] = left_halves[i - 1]*right_halves[n - i]
  end
  return products
end

# Returns the SNF S of G and a map from S to G and one from S to codomain(GtoR).
# If RtoQ is given, then the generators of S are computed modulo codomain(RtoQ).
# It is assumed, that codomain(GtoR) == domain(RtoQ) in this case.
function snf(G::GrpAbFinGen, GtoR::Union{GrpAbFinGenToAbsOrdQuoRingMultMap, GrpAbFinGenToAbsOrdMap}, modulo_map::AbsOrdQuoMap...)
  S, StoG = snf(G)

  R = codomain(GtoR)

  if ngens(G) == 0
    StoR = typeof(GtoR)(S, R, GtoR.generators, GtoR.discrete_logarithm)
    return S, StoG, StoR
  end

  function disclog(x)
    @assert parent(x) === R
    y = GtoR.discrete_logarithm(x)
    a = StoG\(G(y))
    return fmpz[ a[j] for j = 1:ngens(S) ]
  end

  modulo = (length(modulo_map) == 1)
  if modulo
    RtoQ = modulo_map[1]
  end

  gens_snf = Vector{elem_type(R)}(undef, ngens(S))
  for i = 1:ngens(S)
    x = StoG(S[i]).coeff
    for j = 1:ngens(G)
      x[1, j] = mod(x[1, j], S.snf[end])
    end
    modulo ? y = one(codomain(RtoQ)) : y = one(R)
    for j = 1:ngens(G)
      if modulo
        y *= RtoQ(GtoR.generators[j])^(x[1, j])
      else
        y *= GtoR.generators[j]^(x[1, j])
      end
    end
    modulo ? gens_snf[i] = RtoQ\y : gens_snf[i] = y
  end

  if isdefined(GtoR, :modulus)
    StoR = typeof(GtoR)(S, R, gens_snf, disclog, GtoR.modulus)
  else
    StoR = typeof(GtoR)(S, R, gens_snf, disclog)
  end

  if isdefined(GtoR, :tame)
    StoR.tame = GtoR.tame
  end
  if isdefined(GtoR, :wild)
    StoR.wild = GtoR.wild
  end

  return S, StoG, StoR
end

snf(GtoR::Union{GrpAbFinGenToAbsOrdQuoRingMultMap, GrpAbFinGenToAbsOrdMap}, modulo_map::AbsOrdQuoMap...) = snf(domain(GtoR), GtoR, modulo_map...)

# Let QG = codomain(GtoQ) and QH = codomain(HtoQ) with QG = O/I and QH = O/J.
# This function returns the direct product of G and H and a map from the
# product to O/(IJ).
# It is assumed that I and J are coprime.
function direct_product(G::GrpAbFinGen, GtoQ::GrpAbFinGenToAbsOrdQuoRingMultMap, H::GrpAbFinGen, HtoQ::GrpAbFinGenToAbsOrdQuoRingMultMap)
  return direct_product([G, H], [GtoQ, GtoH])
end

# Let G_i be the groups and Q_i be the codomains of the maps such that
# Q_i = O/I_i for ideals I_i.
# This function returns the direct product of the G_i, a map from the
# product to O/(\prod_i I_i) and a map from O to this quotient.
# It is assumed that the ideals I_i are coprime.
function direct_product(groups::Vector{GrpAbFinGen}, maps::Vector{T}) where T <: GrpAbFinGenToAbsOrdQuoRingMultMap
  @assert length(groups) == length(maps)
  @assert length(groups) != 0

  if length(groups) == 1
    Q = codomain(maps[1])
    return groups[1], maps[1], AbsOrdQuoMap(base_ring(Q), Q)
  end

  ideals = Vector{typeof(ideal(codomain(maps[1])))}(undef, length(maps))
  for i = 1:length(maps)
    ideals[i] = ideal(codomain(maps[i]))
  end

  O = order(ideals[1])
  @assert all([ order(ideals[i]) == O for i = 1:length(ideals) ])
  Q, OtoQ = quo(O, prod(ideals))
  return _direct_product(groups, maps, ideals, Q)..., OtoQ
end

# This function returns the direct product of the G_i and a map from the
# product to Q.
function direct_product(groups::Vector{GrpAbFinGen}, maps::Vector{T}, Q::AbsOrdQuoRing) where T <: GrpAbFinGenToAbsOrdQuoRingMultMap
  @assert length(groups) == length(maps)
  @assert length(groups) != 0

  ideals = Vector{ideal_type(base_ring(Q))}(undef, length(maps))
  for i = 1:length(maps)
    ideals[i] = ideal(codomain(maps[i]))
  end
  return _direct_product(groups, maps, ideals, Q)
end

# This function does NOT require length(group) == length(ideals), but only
# length(groups) <= length(ideals). The generators of groups[i] will be made
# coprime to ALL ideals except ideals[i] in any case.
# If tame_wild is true, it is assumed, that for each map in maps the field tame
# is defined and contains exactly one key. The returned map will then contain
# a dictionary of all given tame respectively wild parts with changed generators.
function _direct_product(groups::Vector{GrpAbFinGen}, maps::Vector{U}, ideals::Vector{T}, Q::AbsOrdQuoRing{S, T}, tame_wild::Bool = false) where {S, T, U <: GrpAbFinGenToAbsOrdQuoRingMultMap}
  @assert length(groups) == length(maps)
  @assert length(groups) != 0
  @assert length(ideals) >= length(groups)

  if length(groups) == 1 && length(ideals) == 1
    if codomain(maps[1]) === Q
      return groups[1], maps[1]
    end

    gens = map(Q, [ g.elem for g in maps[1].generators ])

    local disc_log1
    let maps = maps
      function disc_log1(x::AbsOrdQuoRingElem)
        xx = codomain(maps[1])(x.elem)
        return maps[1].discrete_logarithm(xx)
      end
    end

    m = GrpAbFinGenToAbsOrdQuoRingMultMap(groups[1], Q, gens, disc_log1)
    if tame_wild
      m.tame = maps[1].tame
      m.wild = maps[1].wild
    end
    return groups[1], m
  end

  G = direct_product(groups..., task = :none)::GrpAbFinGen

  if tame_wild
    tame = Dict{T, GrpAbFinGenToAbsOrdMap{S}}()
    wild = Dict{T, GrpAbFinGenToAbsOrdMap{S}}()
  end

  O = base_ring(Q)
  oneO = one(O)
  generators = Vector{elem_type(Q)}()
  moduli = _compute_products_for_make_coprime(ideals)
  for i = 1:length(groups)
    if tame_wild
      @assert length(keys(maps[i].tame)) == 1
      p = first(keys(maps[i].tame))
      wild_generators = Vector{elem_type(O)}()
    end
    for j = 1:ngens(groups[i])
      g = crt(maps[i].generators[j].elem, ideals[i], oneO, moduli[i])
      push!(generators, Q(g))
      if tame_wild
        if j == 1
          tame[p] = GrpAbFinGenToAbsOrdMap(domain(maps[i].tame[p]), O, [ g ], maps[i].tame[p].discrete_logarithm)
        else
          push!(wild_generators, g)
        end
      end
    end
    if tame_wild && haskey(maps[i].wild, p)
      wild[p] = GrpAbFinGenToAbsOrdMap(domain(maps[i].wild[p]), O, wild_generators, maps[i].wild[p].discrete_logarithm)
    end
  end

  if length(groups) == 1
    local disc_log2
    let maps = maps
      function disc_log2(x::AbsOrdQuoRingElem)
        xx = codomain(maps[1])(x.elem)
        return maps[1].discrete_logarithm(xx)
      end
    end

    m = GrpAbFinGenToAbsOrdQuoRingMultMap(G, Q, generators, disc_log2)
    if tame_wild
      m.tame = tame
      m.wild = wild
    end
    return G, m
  end

  local disc_log
  let maps = maps
    function disc_log(a::AbsOrdQuoRingElem)
      result = Vector{fmpz}()
      for map in maps
        aa = codomain(map)(a.elem)
        append!(result, map.discrete_logarithm(aa))
      end
      return result
    end
  end

  m = GrpAbFinGenToAbsOrdQuoRingMultMap(G, Q, generators, disc_log)
  if tame_wild
    m.tame = tame
    m.wild = wild
  end
  return G, m
end

function _direct_product!(ideals_and_maps::Vector{Tuple{NfOrdIdl, Vector{GrpAbFinGenToNfAbsOrdMap}}}, Q::NfOrdQuoRing)
  @assert !isempty(ideals_and_maps)

  groups = Vector{GrpAbFinGen}()
  ideals = [x[1] for x in ideals_and_maps]
  O = base_ring(Q)
  oneO = O(1)
  generators = Vector{NfOrdQuoRingElem}()
  if length(ideals) != 1
    moduli = _compute_products_for_make_coprime(ideals)
  end
  for i = 1:length(ideals)
    for map in ideals_and_maps[i][2]
      push!(groups, domain(map))
      for j = 1:length(map.generators)
        if length(ideals) == 1
          push!(generators, Q(map.generators[j]))
        else
          g = crt(map.generators[j], ideals[i], oneO, moduli[i])
          push!(generators, Q(g))
          map.generators[j] = g
        end
      end
    end
  end

  @assert !isempty(groups)
  G = direct_product(groups..., task = :none)

  function disc_log(a::NfOrdQuoRingElem)
    result = Vector{fmpz}()
    for i = 1:length(ideals)
      for map in ideals_and_maps[i][2]
        append!(result, map.discrete_logarithm(a.elem))
      end
    end
    return result
  end

  return G, GrpAbFinGenToAbsOrdQuoRingMultMap(G, Q, generators, disc_log)
end

# This function returns the same as _direct_product called with the values of
# ideals_and_maps as maps and the keys as ideals.
# However, there can be more than one map (or group) for one ideal and the
# generators of the maps are changed IN PLACE.
function _direct_product!(ideals_and_maps::Dict{NfOrdIdl, Vector{GrpAbFinGenToNfAbsOrdMap}}, Q::NfOrdQuoRing)
  @assert !isempty(ideals_and_maps)

  ideals = collect(keys(ideals_and_maps))

  groups = Vector{GrpAbFinGen}()
  O = base_ring(Q)
  oneO = O(1)
  generators = Vector{NfOrdQuoRingElem}()
  if length(ideals) != 1
    moduli = _compute_products_for_make_coprime(ideals)
  end
  for i = 1:length(ideals)
    for map in ideals_and_maps[ideals[i]]
      push!(groups, domain(map))
      for j = 1:length(map.generators)
        if length(ideals) == 1
          push!(generators, Q(map.generators[j]))
        else
          g = crt(map.generators[j], ideals[i], oneO, moduli[i])
          push!(generators, Q(g))
          map.generators[j] = g
        end
      end
    end
  end

  @assert !isempty(groups)
  G = direct_product(groups..., task = :none)

  function disc_log(a::NfOrdQuoRingElem)
    result = Vector{fmpz}()
    for ideal in ideals
      for map in ideals_and_maps[ideal]
        append!(result, map.discrete_logarithm(a.elem))
      end
    end
    return result
  end

  return G, GrpAbFinGenToAbsOrdQuoRingMultMap(G, Q, generators, disc_log)
end


text:
@testset "Mixed Integer Programs: $solver" for solver in solvers
    if can_solve_mip(solver)
        mip_solver = !can_solve_mip(solver) ? GLPKSolverMIP() : solver

        @testset "lp fallback interface" begin
            x = Variable()
            p = minimize(x, x>=4.3)
            @test vexity(p) == AffineVexity()
            solve!(p, solver)
            @test p.optval ≈ 4.3 atol=TOL

            x = Variable(2)
            p = minimize(norm(x,1), x[1]>=4.3)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 4.3 atol=TOL
        end

        @testset "integer variables" begin
            x = Variable(:Int)
            p = minimize(x, x>=4.3)
            @test vexity(p) == AffineVexity()
            solve!(p, mip_solver)
            @test p.optval ≈ 5 atol=TOL

            x = Variable(2, :Int)
            p = minimize(sum(x), x>=4.3)
            @test vexity(p) == AffineVexity()
            solve!(p, mip_solver)
            @test p.optval ≈ 10 atol=TOL

            x = Variable(:Int)
            y = Variable()
            p = minimize(sum(x + y), x>=4.3, y>=7)
            @test vexity(p) == AffineVexity()
            solve!(p, mip_solver)
            @test p.optval ≈ 12 atol=TOL

            x = Variable(2, :Int)
            p = minimize(norm(x, 1), x[1]>=4.3)
            @test vexity(p) == ConvexVexity()
            solve!(p, mip_solver)
            @test p.optval ≈ 5 atol=TOL

            x = Variable(2, :Int)
            p = minimize(sum(x), x[1]>=4.3, x>=0)
            @test vexity(p) == AffineVexity()
            solve!(p, mip_solver)
            @test p.optval ≈ 5 atol=TOL

            x = Variable(2, :Int)
            p = minimize(sum(x), x>=.5)
            @test vexity(p) == AffineVexity()
            solve!(p, mip_solver)
            @test p.optval ≈ 2 atol=TOL
        end

        @testset "binary variables" begin
            x = Variable(2, :Bin)
            p = minimize(sum(x), x>=.5)
            @test vexity(p) == AffineVexity()
            solve!(p, mip_solver)
            @test p.optval ≈ 2 atol=TOL

            x = Variable(2, :Bin)
            p = minimize(sum(x), x[1]>=.5, x>=0)
            @test vexity(p) == AffineVexity()
            solve!(p, mip_solver)
            @test p.optval ≈ 1 atol=TOL
        end
    end
end

ext: jl


text:
export FactorGraph, FactorNode, VarNode, fromUAI
abstract type FGNode end

mutable struct VarNode <: FGNode
    dim:: Int
end

mutable struct FactorNode <: FGNode
    neighbs:: Array{VarNode}
    factor::Dict{Array{Int}, Float64}
end

mutable struct FactorGraph
    facs::Array{FactorNode}
    vars::Array{VarNode}
end

function fromUAI(filename)
    # Load Lines
    lines = collect(eachline(filename))
    # Lines don't actually mean anything, just take whitespace separated
    tokens = filter(x -> length(x) > 0, split(join(lines, " "), r"\s+"))
    numVars = parse(Int, tokens[2])

    # Create variable nodes
    cards = map(x -> parse(Int, x), tokens[3:numVars+2])
    varNodes = [VarNode(c) for c in cards]

    # Create factor nodes
    numFacs = parse(Int, tokens[3+numVars])
    curr_token = 4 + numVars
    # Read in the neighbour list for factors
    neighbs = Array{Array{Int}}(undef, numFacs)
    for fac_ind in 1:numFacs
        num_neighbs = parse(Int, tokens[curr_token])
        # We're 1-indexed in julia
        neighbs[fac_ind] = map(x -> parse(Int, x) + 1, tokens[curr_token+1:curr_token+num_neighbs])
        # Update where we're looking
        curr_token += num_neighbs + 1
    end
    # Read in factor description
    facNodes = Array{FactorNode}(undef, numFacs)
    for fac_ind in 1:numFacs
        # Create the indices for the map
        neighb_cards = map(x -> varNodes[x].dim, neighbs[fac_ind])
        num_fac_entries = reduce(*, neighb_cards)
        # @show neighbs[fac_ind]
        # @show neighb_cards
        labels = reshape(permutedims(collect(Base.product([1:n for n in (neighb_cards)]...)), reverse(1:length(neighb_cards))), num_fac_entries)
        # Convert tuples to arrays
        labels = map(x -> [x...], labels)
        values = map(x -> parse(Float64, x), tokens[curr_token+1:curr_token+num_fac_entries])
        factor = Dict(labels .=> values)
        facNodes[fac_ind] = FactorNode([varNodes[x] for x in neighbs[fac_ind]], factor)
        # Update where we're looking
        curr_token += num_fac_entries + 1
    end

    FactorGraph(facNodes, varNodes)
end


text:
@testset "Stewart Problems" begin
    (a, b, c, n, t, x, y, ) = @variables a b c n t x y

    #= ::Package:: =#

    #= ::Title:: =#
    #=James*Stewart - Calculus*(1987)=#

    #= ::Section::Closed:: =#
    #=Section*7.1 - Integration*by*Parts=#

    @test_int [x^n, x, 1, x^(1 + n)/(1 + n)]
    @test_int [ℯ^x, x, 1, ℯ^x]
    @test_int [1/x, x, 1, log(x)]
    @test_int [a^x, x, 1, a^x/log(a)]

    @test_int [sin(x), x, 1, -cos(x)]
    @test_int [cos(x), x, 1, sin(x)]
    @test_int [sec(x)^2, x, 2, tan(x)]
    @test_int [csc(x)^2, x, 2, -cot(x)]
    @test_int [sec(x)*tan(x), x, 2, sec(x)]
    @test_int [csc(x)*cot(x), x, 2, -csc(x)]

    @test_int [sinh(x), x, 1, cosh(x)]
    @test_int [cosh(x), x, 1, sinh(x)]
    @test_int [tan(x), x, 1, -log(cos(x))]
    @test_int [cot(x), x, 1, log(sin(x))]

    @test_int [x*sin(x), x, 2, -x*cos(x) + sin(x)]
    @test_int [log(x), x, 1, -x + x*log(x)]
    @test_int [x^2*ℯ^x, x, 3, 2*ℯ^x - 2*ℯ^x*x + ℯ^x*x^2]
    @test_int [ℯ^x*sin(x), x, 1, -1/2*ℯ^x*cos(x) + (1/2)*ℯ^x*sin(x)]
    @test_int [atan(x), x, 2, x*atan(x) - log(1 + x^2)/2]

    @test_int [x*ℯ^(2*x), x, 2, -ℯ^(2*x)/4 + (ℯ^(2*x)*x)/2]
    @test_int [x*cos(x), x, 2, cos(x) + x*sin(x)]
    @test_int [x*sin(4*x), x, 2, -(x*cos(4*x))/4 + sin(4*x)/16]
    @test_int [x*log(x), x, 1, -x^2/4 + (x^2*log(x))/2]
    @test_int [x^2*cos(3*x), x, 3, (2*x*cos(3*x))/9 - (2*sin(3*x))/27 + (x^2*sin(3*x))/3]
    @test_int [x^2*sin(2*x), x, 3, cos(2*x)/4 - (x^2*cos(2*x))/2 + (x*sin(2*x))/2]
    @test_int [log(x)^2, x, 2, 2*x - 2*x*log(x) + x*log(x)^2]
    @test_int [asin(x), x, 2, sqrt(1 - x^2) + x*asin(x)]
    @test_int [t*cos(t)*sin(t), t, 3, -t/4 + (cos(t)*sin(t))/4 + (t*sin(t)^2)/2]
    @test_int [t*sec(t)^2, t, 2, log(cos(t)) + t*tan(t)]

    @test_int [t^2*log(t), t, 1, -t^3/9 + (t^3*log(t))/3]
    @test_int [t^3*ℯ^t, t, 4, -6*ℯ^t + 6*ℯ^t*t - 3*ℯ^t*t^2 + ℯ^t*t^3]
    @test_int [ℯ^(2*t)*sin(3*t), t, 1, (-3*ℯ^(2*t)*cos(3*t))/13 + (2*ℯ^(2*t)*sin(3*t))/13]
    @test_int [cos(3*t)/ℯ^t, t, 1, -cos(3*t)/(10*ℯ^t) + (3*sin(3*t))/(10*ℯ^t)]
    @test_int [y*sinh(y), y, 2, y*cosh(y) - sinh(y)]
    @test_int [y*cosh(a*y), y, 2, -(cosh(a*y)/a^2) + (y*sinh(a*y))/a]
    @test_int [t/ℯ^t, t, 2, -ℯ^(-t) - t/ℯ^t]
    @test_int [sqrt(t)*log(t), t, 1, (-4*t^(3/2))/9 + (2*t^(3/2)*log(t))/3]
    @test_int [x*cos(2*x), x, 2, cos(2*x)/4 + (x*sin(2*x))/2]
    @test_int [x^2/ℯ^x, x, 3, -2/ℯ^x - (2*x)/ℯ^x - x^2/ℯ^x]

    @test_int [acos(x), x, 2, -sqrt(1 - x^2) + x*acos(x)]
    @test_int [x*csc(x)^2, x, 2, -(x*cot(x)) + log(sin(x))]
    @test_int [sin(3*x)*cos(5*x), x, 1, (1/4)*cos(2*x) - (1/16)*cos(8*x)]
    @test_int [sin(2*x)*sin(4*x), x, 1, (1/4)*sin(2*x) - (1/12)*sin(6*x)]
    @test_int [cos(x)*log(sin(x)), x, 2, -sin(x) + log(sin(x))*sin(x)]
    @test_int [x^3*ℯ^(x^2), x, 2, -ℯ^x^2/2 + (ℯ^x^2*x^2)/2]
    @test_int [(3 + 2*x)*ℯ^x, x, 2, -2*ℯ^x + ℯ^x*(3 + 2*x)]
    @test_int [x*5^x, x, 2, -(5^x/log(5)^2) + (5^x*x)/log(5)]
    @test_int [cos(log(x)), x, 1, (x*cos(log(x)))/2 + (x*sin(log(x)))/2]
    @test_int [ℯ^sqrt(x), x, 3, -2*ℯ^sqrt(x) + 2*ℯ^sqrt(x)*sqrt(x)]

    @test_int [log(sqrt(x)), x, 1, -x/2 + x*log(sqrt(x))]
    @test_int [sin(log(x)), x, 1, -(x*cos(log(x)))/2 + (x*sin(log(x)))/2]
    @test_int [sin(sqrt(x)), x, 3, -2*sqrt(x)*cos(sqrt(x)) + 2*sin(sqrt(x))]
    @test_int [x^5*cos(x^3), x, 3, cos(x^3)/3 + (x^3*sin(x^3))/3]
    @test_int [x^5*ℯ^(x^2), x, 3, ℯ^x^2 - ℯ^x^2*x^2 + (ℯ^x^2*x^4)/2]
    @test_int [x*atan(x), x, 3, -x/2 + atan(x)/2 + (x^2*atan(x))/2]

    @test_int [x*cos(pi*x), x, 2, cos(pi*x)/pi^2 + (x*sin(pi*x))/pi]
    @test_int [sqrt(x)*log(x), x, 1, (-4*x^(3/2))/9 + (2*x^(3/2)*log(x))/3]

    #= ::Section::Closed:: =#
    #=Section*7.2 - Trigonometric*Integrals=#

    @test_int [sin(3*x)^2, x, 2, x/2 - (cos(3*x)*sin(3*x))/6]
    @test_int [cos(x)^2, x, 2, x/2 + (cos(x)*sin(x))/2]
    @test_int [cos(x)^4, x, 3, (3*x)/8 + (3*cos(x)*sin(x))/8 + (cos(x)^3*sin(x))/4]
    @test_int [sin(x)^3, x, 2, -cos(x) + cos(x)^3/3]
    @test_int [sin(x)^3*cos(x)^4, x, 3, -cos(x)^5/5 + cos(x)^7/7]
    @test_int [sin(x)^4*cos(x)^3, x, 3, sin(x)^5/5 - sin(x)^7/7]
    @test_int [sin(x)^4*cos(x)^2, x, 4, x/16 + (cos(x)*sin(x))/16 - (cos(x)^3*sin(x))/8 - (cos(x)^3*sin(x)^3)/6]
    @test_int [sin(x)^2*cos(x)^2, x, 3, x/8 + (cos(x)*sin(x))/8 - (cos(x)^3*sin(x))/4]
    @test_int [(1 - sin(2*x))^2, x, 1, (3*x)/2 + cos(2*x) - (cos(2*x)*sin(2*x))/4]
    @test_int [sin(x + pi/6)*cos(x), x, 3, x/4 - cos(pi/6 + 2*x)/4]

    @test_int [cos(x)^5*sin(x)^5, x, 4, sin(x)^6/6 - sin(x)^8/4 + sin(x)^10/10]
    @test_int [sin(x)^6, x, 4, (5*x)/16 - (5*cos(x)*sin(x))/16 - (5*cos(x)*sin(x)^3)/24 - (cos(x)*sin(x)^5)/6]
    @test_int [cos(x)^6, x, 4, (5*x)/16 + (5/16)*cos(x)*sin(x) + (5/24)*cos(x)^3*sin(x) + (1/6)*cos(x)^5*sin(x)]
    @test_int [sin(2*x)^2*cos(2*x)^4, x, 4, x/16 + (1/32)*cos(2*x)*sin(2*x) + (1/48)*cos(2*x)^3*sin(2*x) - (1/12)*cos(2*x)^5*sin(2*x)]
    @test_int [sin(x)^5, x, 2, -cos(x) + (2*cos(x)^3)/3 - cos(x)^5/5]
    @test_int [sin(x)^4*cos(x)^4, x, 5, (3*x)/128 + (3/128)*cos(x)*sin(x) + (1/64)*cos(x)^3*sin(x) - (1/16)*cos(x)^5*sin(x) - (1/8)*cos(x)^5*sin(x)^3]
    @test_int [sin(x)^3*sqrt(cos(x)), x, 3, (-2*cos(x)^(3/2))/3 + (2*cos(x)^(7/2))/7]
    @test_int [cos(x)^3*sqrt(sin(x)), x, 3, (2/3)*sin(x)^(3/2) - (2/7)*sin(x)^(7/2)]
    @test_int [cos(sqrt(x))^2/sqrt(x), x, 3, sqrt(x) + cos(sqrt(x))*sin(sqrt(x))]
    @test_int [x*sin(x^2)^3, x, 3, -cos(x^2)/2 + cos(x^2)^3/6]

    @test_int [cos(x)^2*tan(x)^3, x, 3, cos(x)^2/2 - log(cos(x))]
    @test_int [cot(x)^5*sin(x)^2, x, 4, (-(1/2))*csc(x)^2 - 2*log(sin(x)) + sin(x)^2/2]
    @test_int [(1 - sin(x))/cos(x), x, 2, log(1 + sin(x))]
    @test_int [1/(1 - sin(x)), x, 1, cos(x)/(1 - sin(x))]
    @test_int [tan(x)^2, x, 2, -x + tan(x)]
    @test_int [tan(x)^4, x, 3, x - tan(x) + tan(x)^3/3]
    @test_int [sec(x)^4, x, 2, tan(x) + tan(x)^3/3]
    @test_int [sec(x)^6, x, 2, tan(x) + (2*tan(x)^3)/3 + tan(x)^5/5]
    @test_int [tan(x)^4*sec(x)^2, x, 2, tan(x)^5/5]
    @test_int [tan(x)^2*sec(x)^4, x, 3, tan(x)^3/3 + tan(x)^5/5]

    @test_int [tan(x)*sec(x)^3, x, 2, sec(x)^3/3]
    @test_int [tan(x)^3*sec(x)^3, x, 3, -sec(x)^3/3 + sec(x)^5/5]
    @test_int [tan(x)^5, x, 3, -log(cos(x)) - tan(x)^2/2 + tan(x)^4/4]
    @test_int [tan(x)^6, x, 4, -x + tan(x) - tan(x)^3/3 + tan(x)^5/5]
    @test_int [tan(x)^5*sec(x), x, 3, sec(x) - (2*sec(x)^3)/3 + sec(x)^5/5]
    @test_int [tan(x)^5*sec(x)^3, x, 3, sec(x)^3/3 - (2*sec(x)^5)/5 + sec(x)^7/7]
    @test_int [tan(x)*sec(x)^6, x, 2, sec(x)^6/6]
    @test_int [tan(x)^3*sec(x)^6, x, 3, (-(1/6))*sec(x)^6 + sec(x)^8/8]
    @test_int [sec(x)^2/cot(x), x, 2, sec(x)^2/2]
    @test_int [sec(x)*tan(x)^2, x, 2, -atanh(sin(x))/2 + (sec(x)*tan(x))/2]

    @test_int [cot(x)^2, x, 2, -x - cot(x)]
    @test_int [cot(x)^3, x, 2, -cot(x)^2/2 - log(sin(x))]
    @test_int [cot(x)^4*csc(x)^4, x, 3, -cot(x)^5/5 - cot(x)^7/7]
    @test_int [cot(x)^3*csc(x)^4, x, 3, csc(x)^4/4 - csc(x)^6/6]
    @test_int [csc(x), x, 1, -atanh(cos(x))]
    @test_int [csc(x)^3, x, 2, -atanh(cos(x))/2 - (cot(x)*csc(x))/2]
    @test_int [cos(x)^2/sin(x), x, 3, -atanh(cos(x)) + cos(x)]
    @test_int [1/sin(x)^4, x, 2, -cot(x) - cot(x)^3/3]
    @test_int [sin(5*x)*sin(2*x), x, 1, sin(3*x)/6 - sin(7*x)/14]
    @test_int [sin(3*x)*cos(x), x, 1, -cos(2*x)/4 - cos(4*x)/8]

    @test_int [cos(3*x)*cos(4*x), x, 1, sin(x)/2 + sin(7*x)/14]
    @test_int [sin(3*x)*sin(6*x), x, 1, sin(3*x)/6 - sin(9*x)/18]
    @test_int [sin(x)*cos(x)^5, x, 2, (-(1/6))*cos(x)^6]
    @test_int [cos(x)*cos(2*x)*cos(3*x), x, 5, x/4 + (1/8)*sin(2*x) + (1/16)*sin(4*x) + (1/24)*sin(6*x)]
    @test_int [(1 - tan(x)^2)/sec(x)^2, x, 2, cos(x)*sin(x)]
    @test_int [(cos(x) + sin(x))/sin(2*x), x, 6, -atanh(cos(x))/2 + atanh(sin(x))/2]

    @test_int [sin(x)^2*tan(x), x, 3, cos(x)^2/2 - log(cos(x))]
    @test_int [cos(x)^2*cot(x)^3, x, 4, -csc(x)^2/2 - 2*log(sin(x)) + sin(x)^2/2]
    @test_int [sec(x)^3*tan(x), x, 2, sec(x)^3/3]
    @test_int [sec(x)^3*tan(x)^3, x, 3, -sec(x)^3/3 + sec(x)^5/5]

    #= ::Section::Closed:: =#
    #=Section*7.3 - Trigonometric*Substitution=#

    @test_int [sqrt(9 - x^2)/x^2, x, 2, -(sqrt(9 - x^2)/x) - asin(x/3)]
    @test_int [1/(x^2*sqrt(4 + x^2)), x, 1, -sqrt(4 + x^2)/(4*x)]
    @test_int [x/sqrt(4 + x^2), x, 1, sqrt(4 + x^2)]
    @test_int [1/sqrt(-a^2 + x^2), x, 2, atanh(x/sqrt(-a^2 + x^2))]
    @test_int [x^3/(9 + 4*x^2)^(3/2), x, 3, 9/(16*sqrt(9 + 4*x^2)) + sqrt(9 + 4*x^2)/16]
    @test_int [x/sqrt(3 - 2*x - x^2), x, 3, -sqrt(3 - 2*x - x^2) + asin((-1 - x)/2)]

    @test_int [1/(x^2*sqrt(1 - x^2)), x, 1, -(sqrt(1 - x^2)/x)]
    @test_int [x^3*sqrt(4 - x^2), x, 3, (-4*(4 - x^2)^(3/2))/3 + (4 - x^2)^(5/2)/5]
    @test_int [x/sqrt(1 - x^2), x, 1, -sqrt(1 - x^2)]
    @test_int [x*sqrt(4 - x^2), x, 1, -(4 - x^2)^(3/2)/3]
    @test_int [sqrt(1 - 4*x^2), x, 2, (x*sqrt(1 - 4*x^2))/2 + asin(2*x)/4]
    @test_int [x^3/sqrt(x^2 + 4), x, 3, -4*sqrt(4 + x^2) + (1/3)*(4 + x^2)^(3/2)]
    @test_int [1/sqrt(9 + x^2), x, 1, asinh(x/3)]
    @test_int [sqrt(x^2 + 1), x, 2, (x*sqrt(1 + x^2))/2 + asinh(x)/2]
    @test_int [1/(x^3*sqrt(x^2 - 16)), x, 4, sqrt(-16 + x^2)/(32*x^2) + atan(sqrt(-16 + x^2)/4)/128]
    @test_int [sqrt(x^2 - a^2)/x^4, x, 1, (-a^2 + x^2)^(3/2)/(3*a^2*x^3)]

    @test_int [sqrt(9*x^2 - 4)/x, x, 4, sqrt(-4 + 9*x^2) - 2*atan(sqrt(-4 + 9*x^2)/2)]
    @test_int [1/(x^2*sqrt(16*x^2 - 9)), x, 1, sqrt(-9 + 16*x^2)/(9*x)]
    @test_int [x^2/(a^2 - x^2)^(3/2), x, 3, x/sqrt(a^2 - x^2) - atan(x/sqrt(a^2 - x^2))]
    @test_int [x^2/sqrt(5 - x^2), x, 2, -(x*sqrt(5 - x^2))/2 + (5*asin(x/sqrt(5)))/2]
    @test_int [1/(x*sqrt(3 + x^2)), x, 3, -(atanh(sqrt(3 + x^2)/sqrt(3))/sqrt(3))]
    @test_int [x/(x^2 + 4)^(5/2), x, 1, -1/(3*(4 + x^2)^(3/2))]
    @test_int [x^3*sqrt(4 - 9*x^2), x, 3, (-4*(4 - 9*x^2)^(3/2))/243 + (4 - 9*x^2)^(5/2)/405]
    @test_int [x^2*sqrt(9 - x^2), x, 3, (-9*x*sqrt(9 - x^2))/8 + (x^3*sqrt(9 - x^2))/4 + (81*asin(x/3))/8]
    @test_int [5*x*sqrt(1 + x^2), x, 2, (5*(1 + x^2)^(3/2))/3]
    @test_int [1/(4*x^2 - 25)^(3/2), x, 1, -x/(25*sqrt(-25 + 4*x^2))]

    @test_int [sqrt(2*x - x^2), x, 3, (-(1/2))*(1 - x)*sqrt(2*x - x^2) - (1/2)*asin(1 - x)]
    @test_int [1/sqrt(x^2 + 4*x + 8), x, 2, asinh((2 + x)/2)]
    @test_int [1/sqrt(9*x^2 + 6*x - 8), x, 2, atanh((1 + 3*x)/sqrt(-8 + 6*x + 9*x^2))/3]
    @test_int [x^2/sqrt(4*x - x^2), x, 4, -3*sqrt(4*x - x^2) - (1/2)*x*sqrt(4*x - x^2) - 6*asin(1 - x/2)]
    @test_int [1/(2 + 2*x + x^2)^2, x, 3, (1 + x)/(2*(2 + 2*x + x^2)) + atan(1 + x)/2]
    @test_int [1/(5 - 4*x - x^2)^(5/2), x, 2, (2 + x)/(27*(5 - 4*x - x^2)^(3/2)) + (2*(2 + x))/(243*sqrt(5 - 4*x - x^2))]
    @test_int [ℯ^t*sqrt(9 - ℯ^(2*t)), t, 3, (ℯ^t*sqrt(9 - ℯ^(2*t)))/2 + (9*asin(ℯ^t/3))/2]
    @test_int [sqrt(ℯ^(2*t) - 9), t, 4, sqrt(-9 + ℯ^(2*t)) - 3*atan(sqrt(-9 + ℯ^(2*t))/3)]

    #= ::Section::Closed:: =#
    #=Section*7.4 - Integration*of*Rational*Functions*by*Partial*Fractions=#

    @test_int [1/sqrt(a^2 + x^2), x, 2, atanh(x/sqrt(a^2 + x^2))]
    @test_int [(5 + x)/(-2 + x + x^2), x, 3, 2*log(1 - x) - log(2 + x)]
    @test_int [(x + x^3)/(-1 + x), x, 3, 2*x + x^2/2 + x^3/3 + 2*log(1 - x)]
    @test_int [(-1 + 2*x + x^2)/(-2*x + 3*x^2 + 2*x^3), x, 3, log(1 - 2*x)/10 + log(x)/2 - log(2 + x)/10]
    @test_int [(1 + 4*x - 2*x^2 + x^4)/(1 - x - x^2 + x^3), x, 2, 2/(1 - x) + x + x^2/2 + log(1 - x) - log(1 + x)]
    @test_int [(4 - x + 2*x^2)/(4*x + x^3), x, 6, -atan(x/2)/2 + log(x) + log(4 + x^2)/2]
    @test_int [(2 - 3*x + 4*x^2)/(3 - 4*x + 4*x^2), x, 6, x + atan((1 - 2*x)/sqrt(2))/(4*sqrt(2)) + log(3 - 4*x + 4*x^2)/8]
    @test_int [(1 + x^2 + x^3)/((-1 + x)*x*(1 + x^2)^3*(1 + x + x^2)), x, 14, (1 + x)/(8*(1 + x^2)^2) - (3*(1 - x))/(8*(1 + x^2)) + (3*x)/(16*(1 + x^2)) + (7*atan(x))/16 - atan((1 + 2*x)/sqrt(3))/sqrt(3) + log(1 - x)/8 - log(x) + (15*log(1 + x^2))/16 - log(1 + x + x^2)/2]
    @test_int [(1 - 3*x + 2*x^2 - x^3)/(x*(x^2 + 1)^2), x, 6, -((1 + 2*x)/(2*(1 + x^2))) - 2*atan(x) + log(x) - (1/2)*log(1 + x^2)]
    @test_int [1/(x^2 + 1)^2, x, 2, x/(2*(1 + x^2)) + atan(x)/2]

    @test_int [1/((x - 1)*(2 + x)), x, 3, log(1 - x)/3 - log(2 + x)/3]
    @test_int [7/(-12 + 5*x + 2*x^2), x, 4, (7*log(3 - 2*x))/11 - (7*log(4 + x))/11]
    @test_int [(-4 + 3*x + x^2)/((-1 + 2*x)^2*(3 + 2*x)), x, 2, -9/(32*(1 - 2*x)) + (41*log(1 - 2*x))/128 - (25*log(3 + 2*x))/128]
    @test_int [(-x^2 + x^3)/((-6 + x)*(3 + 5*x)^3), x, 3, -12/(1375*(3 + 5*x)^2) + 201/(15125*(3 + 5*x)) + (20*log(6 - x))/3993 + (1493*log(3 + 5*x))/499125]
    @test_int [1/(-x^3 + x^4), x, 3, 1/(2*x^2) + 1/x + log(1 - x) - log(x)]
    @test_int [(1 - x - x^2 + x^3 + x^4)/(-x + x^3), x, 4, x + x^2/2 - log(x) + log(1 - x^2)/2]

    @test_int [(x^2 - 2)/(x*(x^2 + 2)), x, 3, -log(x) + log(2 + x^2)]
    @test_int [(2 - 4*x^2 + x^3)/((1 + x^2)*(2 + x^2)), x, 8, 6*atan(x) - 5*sqrt(2)*atan(x/sqrt(2)) - log(1 + x^2)/2 + log(2 + x^2)]
    @test_int [(1 + x^2 + x^4)/((1 + x^2)*(4 + x^2)^2), x, 6, (-13*x)/(24*(4 + x^2)) + (25*atan(x/2))/144 + atan(x)/9]
    @test_int [(1 + 16*x)/((5 + x)^2*(-3 + 2*x)*(1 + x + x^2)), x, 6, -79/(273*(5 + x)) + (451*atan((1 + 2*x)/sqrt(3)))/(2793*sqrt(3)) + (200*log(3 - 2*x))/3211 + (2731*log(5 + x))/24843 - (481*log(1 + x + x^2))/5586]
    @test_int [x^4/(9 + x^2)^3, x, 3, -x^3/(4*(9 + x^2)^2) - (3*x)/(8*(9 + x^2)) + atan(x/3)/8]
    @test_int [(19*x)/((-1 + x)^3*(3 + 5*x + 4*x^2)^2), x, 8, -399/(736*(1 - x)^2) - 1843/(4416*(1 - x)) + (19*(39 + 44*x))/(276*(1 - x)^2*(3 + 5*x + 4*x^2)) + (114437*atan((5 + 8*x)/sqrt(23)))/(52992*sqrt(23)) + (209*log(1 - x))/2304 - (209*log(3 + 5*x + 4*x^2))/4608]
    @test_int [(1 + x^2 + x^3)/(2*x^2 + x^3 + x^4), x, 7, -1/(2*x) + atan((1 + 2*x)/sqrt(7))/(4*sqrt(7)) - log(x)/4 + (5*log(2 + x + x^2))/8]
    @test_int [1/(-x^3 + x^6), x, 8, 1/(2*x^2) - atan((1 + 2*x)/sqrt(3))/sqrt(3) + log(1 - x)/3 - log(1 + x + x^2)/6]
    @test_int [x^2/(1 + x), x, 2, -x + x^2/2 + log(1 + x)]
    @test_int [x/(-5 + x), x, 2, x + 5*log(5 - x)]

    @test_int [(-1 + 4*x)/((-1 + x)*(2 + x)), x, 2, log(1 - x) + 3*log(2 + x)]
    @test_int [1/((1 + x)*(2 + x)), x, 3, log(1 + x) - log(2 + x)]
    @test_int [(-5 + 6*x)/(3 + 2*x), x, 2, 3*x - 7*log(3 + 2*x)]
    @test_int [1/((a + x)*(b + x)), x, 3, -(log(a + x)/(a - b)) + log(b + x)/(a - b)]
    @test_int [(1 + x^2)/(-x + x^2), x, 3, x + 2*log(1 - x) - log(x)]
    @test_int [(1 - 12*x + x^2 + x^3)/(-12 + x + x^2), x, 5, x^2/2 + log(3 - x)/7 - log(4 + x)/7]
    @test_int [(3 + 2*x)/(1 + x)^2, x, 2, -(1 + x)^(-1) + 2*log(1 + x)]
    @test_int [1/(x*(1 + x)*(3 + 2*x)), x, 2, log(x)/3 - log(1 + x) + (2*log(3 + 2*x))/3]

    @test_int [(-3 + 5*x + 6*x^2)/(-3*x + 2*x^2 + x^3), x, 3, 2*log(1 - x) + log(x) + 3*log(3 + x)]
    @test_int [x/(4 + 4*x + x^2), x, 3, 2/(2 + x) + log(2 + x)]
    @test_int [1/((-1 + x)^2*(4 + x)), x, 2, 1/(5*(1 - x)) - log(1 - x)/25 + log(4 + x)/25]
    @test_int [x^2/((-3 + x)*(2 + x)^2), x, 2, 4/(5*(2 + x)) + (9*log(3 - x))/25 + (16*log(2 + x))/25]
    @test_int [(-2 + 3*x + 5*x^2)/(2*x^2 + x^3), x, 3, x^(-1) + 2*log(x) + 3*log(2 + x)]
    @test_int [(18 - 2*x - 4*x^2)/(-6 + x + 4*x^2 + x^3), x, 2, log(1 - x) - 2*log(2 + x) - 3*log(3 + x)]
    @test_int [(2*x + x^2)/(4 + 3*x^2 + x^3), x, 1, log(4 + 3*x^2 + x^3)/3]
    @test_int [1/((-1 + x)^2*x^2), x, 2, (1 - x)^(-1) - x^(-1) - 2*log(1 - x) + 2*log(x)]

    @test_int [x^2/(1 + x)^3, x, 2, -1/(2*(1 + x)^2) + 2/(1 + x) + log(1 + x)]
    @test_int [1/(-x^2 + x^4), x, 3, x^(-1) - atanh(x)]
    @test_int [(-x + 2*x^3)/(1 - x^2 + x^4), x, 1, log(1 - x^2 + x^4)/2]
    @test_int [x^3/(1 + x^2), x, 3, x^2/2 - log(1 + x^2)/2]
    @test_int [(-1 + x)/(2 + 2*x + x^2), x, 4, -2*atan(1 + x) + log(2 + 2*x + x^2)/2]

    @test_int [x/(1 + x + x^2), x, 4, -(atan((1 + 2*x)/sqrt(3))/sqrt(3)) + log(1 + x + x^2)/2]
    @test_int [(7 + 5*x + 4*x^2)/(5 + 4*x + 4*x^2), x, 6, x + (3/8)*atan(1/2 + x) + (1/8)*log(5 + 4*x + 4*x^2)]
    @test_int [(5 - 4*x + 3*x^2)/((-1 + x)*(1 + x^2)), x, 5, -3*atan(x) + 2*log(1 - x) + log(1 + x^2)/2]
    @test_int [(3 + 2*x)/(3*x + x^3), x, 6, (2*atan(x/sqrt(3)))/sqrt(3) + log(x) - log(3 + x^2)/2]
    @test_int [1/(-1 + x^3), x, 6, -(atan((1 + 2*x)/sqrt(3))/sqrt(3)) + log(1 - x)/3 - log(1 + x + x^2)/6]
    @test_int [x^3/(1 + x^3), x, 7, x + atan((1 - 2*x)/sqrt(3))/sqrt(3) - log(1 + x)/3 + log(1 - x + x^2)/6]
    @test_int [(-1 - 2*x + x^2)/((-1 + x)^2*(1 + x^2)), x, 5, (-1 + x)^(-1) + atan(x) + log(1 - x) - log(1 + x^2)/2]
    @test_int [x^4/(-1 + x^4), x, 4, x - atan(x)/2 - atanh(x)/2]

    @test_int [(-4 + 6*x - x^2 + 3*x^3)/((1 + x^2)*(2 + x^2)), x, 6, -3*atan(x) + sqrt(2)*atan(x/sqrt(2)) + (3*log(1 + x^2))/2]
    @test_int [(1 + x - 2*x^2 + x^3)/(4 + 5*x^2 + x^4), x, 7, (-3*atan(x/2))/2 + atan(x) + log(4 + x^2)/2]
    @test_int [(-3 + x)/(4 + 2*x + x^2)^2, x, 3, -(7 + 4*x)/(6*(4 + 2*x + x^2)) - (2*atan((1 + x)/sqrt(3)))/(3*sqrt(3))]
    @test_int [(1 + x^4)/(x*(1 + x^2)^2), x, 3, (1 + x^2)^(-1) + log(x)]
    @test_int [(cos(x)*(-3 + 2*sin(x)))/(2 - 3*sin(x) + sin(x)^2), x, 2, log(2 - 3*sin(x) + sin(x)^2)]
    @test_int [(cos(x)^2*sin(x))/(5 + cos(x)^2), x, 3, sqrt(5)*atan(cos(x)/sqrt(5)) - cos(x)]

    @test_int [1/(x^2 + 2*x - 3), x, 3, log(1 - x)/4 - log(3 + x)/4]
    @test_int [1/(x^2 - 2*x), x, 1, log(2 - x)/2 - log(x)/2]
    @test_int [(2*x + 1)/(4*x^2 + 12*x - 7), x, 3, log(1 - 2*x)/8 + (3*log(7 + 2*x))/8]
    @test_int [x/(x^2 + x - 1), x, 3, ((5 - sqrt(5))*log(1 - sqrt(5) + 2*x))/10 + ((5 + sqrt(5))*log(1 + sqrt(5) + 2*x))/10]

    #= ::Section::Closed:: =#
    #=Section*7.5 - Rationalization*Substitutions=#

    @test_int [(-32 + 5*x - 27*x^2 + 4*x^3)/(-70 - 299*x - 286*x^2 + 50*x^3 - 13*x^4 + 30*x^5), x, 6, (3988*atan((1 + 2*x)/sqrt(19)))/(13685*sqrt(19)) - (3146*log(7 - 3*x))/80155 - (334*log(1 + 2*x))/323 + (4822*log(2 + 5*x))/4879 + (11049*log(5 + x + x^2))/260015]
    @test_int [(8 - 13*x^2 - 7*x^3 + 12*x^5)/(4 - 20*x + 41*x^2 - 80*x^3 + 116*x^4 - 80*x^5 + 100*x^6), x, 7, 5828/(9075*(2 - 5*x)) - (313 + 502*x)/(1452*(1 + 2*x^2)) - (251*atan(sqrt(2)*x))/(726*sqrt(2)) + (272*sqrt(2)*atan(sqrt(2)*x))/1331 - (59096*log(2 - 5*x))/99825 + (2843*log(1 + 2*x^2))/7986]
    @test_int [sqrt(4 + x)/x, x, 3, 2*sqrt(4 + x) - 4*atanh(sqrt(4 + x)/2)]
    @test_int [1/(-x^(-1/3) + sqrt(x)), x, 9, 2*sqrt(x) + (3/5)*sqrt(2*(5 - sqrt(5)))*atan((1 - sqrt(5) + 4*x^(1/6))/sqrt(2*(5 + sqrt(5)))) - (3/5)*sqrt(2*(5 + sqrt(5)))*atan((1/2)*sqrt((1/10)*(5 + sqrt(5)))*(1 + sqrt(5) + 4*x^(1/6))) + (6/5)*log(1 - x^(1/6)) - (3/10)*(1 + sqrt(5))*log(2 + x^(1/6) - sqrt(5)*x^(1/6) + 2*x^(1/3)) - (3/10)*(1 - sqrt(5))*log(2 + x^(1/6) + sqrt(5)*x^(1/6) + 2*x^(1/3))]
    @test_int [1/(-4*cos(x) + 3*sin(x)), x, 2, -atanh((3*cos(x) + 4*sin(x))/5)/5]

    @test_int [1/(1 + sqrt(x)), x, 3, 2*sqrt(x) - 2*log(1 + sqrt(x))]
    @test_int [1/(1 + x^(-1/3)), x, 3, 3*x^(1/3) - (3*x^(2/3))/2 + x - 3*log(1 + 1/x^(1/3)) - log(x)]
    @test_int [sqrt(x)/(1 + x), x, 3, 2*sqrt(x) - 2*atan(sqrt(x))]
    @test_int [1/(x*sqrt(1 + x)), x, 2, -2*atanh(sqrt(1 + x))]
    @test_int [1/(-x^(1/3) + x), x, 2, (3*log(1 - x^(2/3)))/2]
    @test_int [1/(x - sqrt(2 + x)), x, 4, (4*log(2 - sqrt(2 + x)))/3 + (2*log(1 + sqrt(2 + x)))/3]
    @test_int [x^2/sqrt(-1 + x), x, 2, 2*sqrt(-1 + x) + (4*(-1 + x)^(3/2))/3 + (2*(-1 + x)^(5/2))/5]
    @test_int [sqrt(-1 + x)/(1 + x), x, 3, 2*sqrt(-1 + x) - 2*sqrt(2)*atan(sqrt(-1 + x)/sqrt(2))]
    @test_int [1/sqrt(1 + sqrt(x)), x, 3, -4*sqrt(1 + sqrt(x)) + (4*(1 + sqrt(x))^(3/2))/3]
    @test_int [sqrt(x)/(x + x^2), x, 3, 2*atan(sqrt(x))]
    @test_int [(1 + sqrt(x))/(-1 + sqrt(x)), x, 3, 4*sqrt(x) + x + 4*log(1 - sqrt(x))]
    @test_int [(1 + x^(-1/3))/(-1 + x^(-1/3)), x, 4, -6*x^(1/3) - 3*x^(2/3) - x - 6*log(1 - x^(1/3))]
    @test_int [x^3/(1 + x^2)^(1/3), x, 3, (-(3/4))*(1 + x^2)^(2/3) + (3/10)*(1 + x^2)^(5/3)]
    @test_int [sqrt(x)/(-x^(-1/3) + sqrt(x)), x, 10, 6*x^(1/6) + x - (3/5)*sqrt(2*(5 + sqrt(5)))*atan((1 - sqrt(5) + 4*x^(1/6))/sqrt(2*(5 + sqrt(5)))) - (3/5)*sqrt(2*(5 - sqrt(5)))*atan((1/2)*sqrt((1/10)*(5 + sqrt(5)))*(1 + sqrt(5) + 4*x^(1/6))) + (6/5)*log(1 - x^(1/6)) - (3/10)*(1 - sqrt(5))*log(2 + x^(1/6) - sqrt(5)*x^(1/6) + 2*x^(1/3)) - (3/10)*(1 + sqrt(5))*log(2 + x^(1/6) + sqrt(5)*x^(1/6) + 2*x^(1/3))]
    @test_int [1/(x^(-1/4) + sqrt(x)), x, 9, 2*sqrt(x) + (4*atan((1 - 2*x^(1/4))/sqrt(3)))/sqrt(3) + (4*log(1 + x^(1/4)))/3 - (2*log(1 - x^(1/4) + sqrt(x)))/3]
    @test_int [1/(x^(-1/3) + x^(-1/4)), x, 4, 12*x^(1/12) - 6*x^(1/6) + 4*x^(1/4) - 3*x^(1/3) + (12*x^(5/12))/5 - 2*sqrt(x) + (12*x^(7/12))/7 - (3*x^(2/3))/2 + (4*x^(3/4))/3 - (6*x^(5/6))/5 + (12*x^(11/12))/11 - x + (12*x^(13/12))/13 - (6*x^(7/6))/7 + (4*x^(5/4))/5 - 12*log(1 + x^(1/12))]
    @test_int [sqrt((1 - x)/x), x, 5, sqrt(-1 + x^(-1))*x - atan(sqrt(-1 + x^(-1)))]
    @test_int [cos(x)/(sin(x) + sin(x)^2), x, 2, log(sin(x)) - log(1 + sin(x))]
    @test_int [ℯ^(2*x)/(2 + 3*ℯ^x + ℯ^(2*x)), x, 4, -log(1 + ℯ^x) + 2*log(2 + ℯ^x)]
    @test_int [1/sqrt(1 + ℯ^x), x, 3, -2*atanh(sqrt(1 + ℯ^x))]
    @test_int [sqrt(1 - ℯ^x), x, 4, 2*sqrt(1 - ℯ^x) - 2*atanh(sqrt(1 - ℯ^x))]
    @test_int [1/(3 - 5*sin(x)), x, 4, -log(cos(x/2) - 3*sin(x/2))/4 + log(3*cos(x/2) - sin(x/2))/4]
    @test_int [1/(cos(x) + sin(x)), x, 2, -(atanh((cos(x) - sin(x))/sqrt(2))/sqrt(2))]
    @test_int [1/(1 - cos(x) + sin(x)), x, 2, -log(1 + cot(x/2))]
    @test_int [1/(4*cos(x) + 3*sin(x)), x, 2, -atanh((3*cos(x) - 4*sin(x))/5)/5]
    @test_int [1/(sin(x) + tan(x)), x, 6, -atanh(cos(x))/2 + (cot(x)*csc(x))/2 - csc(x)^2/2]
    @test_int [1/(2*sin(x) + sin(2*x)), x, 4, log(tan(x/2))/4 + tan(x/2)^2/8]
    @test_int [sec(x)/(1 + sin(x)), x, 4, atanh(sin(x))/2 - 1/(2*(1 + sin(x)))]
    @test_int [1/(b*cos(x) + a*sin(x)), x, 2, -(atanh((a*cos(x) - b*sin(x))/sqrt(a^2 + b^2))/sqrt(a^2 + b^2))]
    @test_int [1/(b^2*cos(x)^2 + a^2*sin(x)^2), x, 2, atan((a*tan(x))/b)/(a*b)]

    #= ::Section::Closed:: =#
    #=Section*7.6 - Strategy*for*Integration=#

    @test_int [x/(-1 + x^2), x, 1, log(1 - x^2)/2]
    @test_int [(1 + sqrt(x))*sqrt(x), x, 2, (2*x^(3/2))/3 + x^2/2]
    @test_int [1/(1 - cos(x)), x, 1, -(sin(x)/(1 - cos(x)))]
    @test_int [sec(x)*tan(x)^2, x, 2, -atanh(sin(x))/2 + (sec(x)*tan(x))/2]
    @test_int [sec(x)^3*tan(x)^3, x, 3, -sec(x)^3/3 + sec(x)^5/5]
    @test_int [ℯ^sqrt(x), x, 3, -2*ℯ^sqrt(x) + 2*ℯ^sqrt(x)*sqrt(x)]
    @test_int [(1 + x^5)/(-10*x - 3*x^2 + x^3), x, 3, 19*x + (3*x^2)/2 + x^3/3 + (3126*log(5 - x))/35 - log(x)/10 - (31*log(2 + x))/14]
    @test_int [1/(x*sqrt(log(x))), x, 2, 2*sqrt(log(x))]
    @test_int [(5 + 2*x)/(-3 + x), x, 2, 2*x + 11*log(3 - x)]
    @test_int [ℯ^(ℯ^x + x), x, 2, ℯ^ℯ^x]

    @test_int [cos(x)^2*sin(x)^2, x, 3, x/8 + (cos(x)*sin(x))/8 - (cos(x)^3*sin(x))/4]
    @test_int [(-cos(x) + sin(x))/(cos(x) + sin(x)), x, 1, -log(cos(x) + sin(x))]
    @test_int [x/sqrt(1 - x^2), x, 1, -sqrt(1 - x^2)]
    @test_int [x^3*log(x), x, 1, -x^4/16 + (x^4*log(x))/4]
    @test_int [sqrt(-2 + x)/(2 + x), x, 3, 2*sqrt(-2 + x) - 4*atan(sqrt(-2 + x)/2)]
    @test_int [x/(2 + x)^2, x, 2, 2/(2 + x) + log(2 + x)]
    @test_int [log(1 + x^2), x, 3, -2*x + 2*atan(x) + x*log(1 + x^2)]
    @test_int [sqrt(1 + log(x))/(x*log(x)), x, 4, -2*atanh(sqrt(1 + log(x))) + 2*sqrt(1 + log(x))]
    @test_int [(1 + sqrt(x))^8, x, 3, (-2*(1 + sqrt(x))^9)/9 + (1 + sqrt(x))^10/5]
    @test_int [sec(x)^4*tan(x)^3, x, 3, (-(1/4))*sec(x)^4 + sec(x)^6/6]

    @test_int [x/(2 - 2*x + x^2), x, 4, -atan(1 - x) + log(2 - 2*x + x^2)/2]
    @test_int [x*asin(x), x, 3, (x*sqrt(1 - x^2))/4 - asin(x)/4 + (x^2*asin(x))/2]
    @test_int [sqrt(9 - x^2)/x, x, 4, sqrt(9 - x^2) - 3*atanh(sqrt(9 - x^2)/3)]
    @test_int [x/(2 + 3*x + x^2), x, 3, -log(1 + x) + 2*log(2 + x)]
    @test_int [x^2*cosh(x), x, 3, -2*x*cosh(x) + 2*sinh(x) + x^2*sinh(x)]
    @test_int [(1 + x + x^3)/(4*x + 2*x^2 + x^4), x, 1, log(4*x + 2*x^2 + x^4)/4]
    @test_int [cos(x)/(1 + sin(x)^2), x, 2, atan(sin(x))]
    @test_int [cos(sqrt(x)), x, 3, 2*cos(sqrt(x)) + 2*sqrt(x)*sin(sqrt(x))]
    @test_int [sin(pi*x), x, 1, -(cos(pi*x)/pi)]
    @test_int [ℯ^(2*x)/(1 + ℯ^x), x, 3, ℯ^x - log(1 + ℯ^x)]

    @test_int [ℯ^(3*x)*cos(5*x), x, 1, (3*ℯ^(3*x)*cos(5*x))/34 + (5*ℯ^(3*x)*sin(5*x))/34]
    @test_int [cos(3*x)*cos(5*x), x, 1, sin(2*x)/4 + sin(8*x)/16]
    @test_int [1/(1 + x + x^2 + x^3), x, 5, atan(x)/2 + log(1 + x)/2 - log(1 + x^2)/4]
    @test_int [x^2*log(1 + x), x, 3, -x/3 + x^2/6 - x^3/9 + log(1 + x)/3 + (x^3*log(1 + x))/3]
    @test_int [x^5/ℯ^x^3, x, 2, -1/(3*ℯ^x^3) - x^3/(3*ℯ^x^3)]
    @test_int [tan(4*x)^2, x, 2, -x + tan(4*x)/4]
    @test_int [1/sqrt(-5 + 12*x + 9*x^2), x, 2, atanh((2 + 3*x)/sqrt(-5 + 12*x + 9*x^2))/3]
    @test_int [x^2*atan(x), x, 4, -x^2/6 + (x^3*atan(x))/3 + log(1 + x^2)/6]
    @test_int [(1 - sqrt(x))/x^(1/3), x, 2, (3*x^(2/3))/2 - (6*x^(7/6))/7]
    @test_int [1/(-ℯ^(-x) + ℯ^x), x, 2, -atanh(ℯ^x)]
    @test_int [x/(10 + 2*x^2 + x^4), x, 3, (1/6)*atan((1/3)*(1 + x^2))]

    @test_int [1/(x^(-1/3) + x), x, 2, (3*log(1 + x^(4/3)))/4]
    @test_int [cos(x)^4*sin(x)^2, x, 4, x/16 + (cos(x)*sin(x))/16 + (cos(x)^3*sin(x))/24 - (cos(x)^5*sin(x))/6]
    @test_int [1/sqrt(5 - 4*x - x^2), x, 2, -asin((-2 - x)/3)]
    @test_int [x/(1 - x^2 + sqrt(1 - x^2)), x, 3, -log(1 + sqrt(1 - x^2))]
    @test_int [(1 + cos(x))*csc(x), x, 2, log(1 - cos(x))]
    @test_int [ℯ^x/(-1 + ℯ^(2*x)), x, 2, -atanh(ℯ^x)]
    @test_int [1/(-8 + x^3), x, 6, -atan((1 + x)/sqrt(3))/(4*sqrt(3)) + log(2 - x)/12 - log(4 + 2*x + x^2)/24]
    @test_int [x^5*cosh(x), x, 6, -120*cosh(x) - 60*x^2*cosh(x) - 5*x^4*cosh(x) + 120*x*sinh(x) + 20*x^3*sinh(x) + x^5*sinh(x)]
    @test_int [log(tan(x))/(sin(x)*cos(x)), x, 1, log(tan(x))^2/2]

    @test_int [-2*x + x^2 + x^3, x, 1, -x^2 + x^3/3 + x^4/4]
    @test_int [(1 + ℯ^x)/(1 - ℯ^x), x, 3, x - 2*log(1 - ℯ^x)]
    @test_int [x/((1 + x^2)*(4 + x^2)), x, 4, log(1 + x^2)/6 - log(4 + x^2)/6]
    @test_int [1/(4 - 5*sin(x)), x, 4, -log(cos(x/2) - 2*sin(x/2))/3 + log(2*cos(x/2) - sin(x/2))/3]
    @test_int [x*(c + x)^(1/3), x, 2, (-3*c*(c + x)^(4/3))/4 + (3*(c + x)^(7/3))/7]
    @test_int [ℯ^x^(1/3), x, 4, 6*ℯ^x^(1/3) - 6*ℯ^x^(1/3)*x^(1/3) + 3*ℯ^x^(1/3)*x^(2/3)]
    @test_int [1/(4 + x + sqrt(1 + x)), x, 5, (-2*atan((1 + 2*sqrt(1 + x))/sqrt(11)))/sqrt(11) + log(4 + x + sqrt(1 + x))]
    @test_int [(1 + x^3)/(-x^2 + x^3), x, 3, 1/x + x + 2*log(1 - x) - log(x)]

    @test_int [(-3 + 4*x + x^2)*sin(2*x), x, 8, (7*cos(2*x))/4 - 2*x*cos(2*x) - (x^2*cos(2*x))/2 + sin(2*x) + (x*sin(2*x))/2]
    @test_int [cos(cos(x))*sin(x), x, 2, -sin(cos(x))]
    @test_int [1/sqrt(16 - x^2), x, 1, asin(x/4)]
    @test_int [x^3/(1 + x)^10, x, 2, 1/(9*(1 + x)^9) - 3/(8*(1 + x)^8) + 3/(7*(1 + x)^7) - 1/(6*(1 + x)^6)]
    @test_int [cot(2*x)^3*csc(2*x)^3, x, 3, csc(2*x)^3/6 - csc(2*x)^5/10]
    @test_int [(x + sin(x))^2, x, 6, x/2 + x^3/3 - 2*x*cos(x) + 2*sin(x) - (cos(x)*sin(x))/2]
    @test_int [ℯ^atan(x)/(1 + x^2), x, 1, ℯ^atan(x)]
    @test_int [1/(x*(1 + x^4)), x, 4, log(x) - log(1 + x^4)/4]
    @test_int [t^3/ℯ^(2*t), t, 4, -3/(8*ℯ^(2*t)) - (3*t)/(4*ℯ^(2*t)) - (3*t^2)/(4*ℯ^(2*t)) - t^3/(2*ℯ^(2*t))]
    @test_int [sqrt(t)/(1 + t^(1/3)), t, 7, -6*t^(1/6) + 2*sqrt(t) - (6*t^(5/6))/5 + (6*t^(7/6))/7 + 6*atan(t^(1/6))]

    @test_int [sin(x)*sin(2*x)*sin(3*x), x, 5, -cos(2*x)/8 - cos(4*x)/16 + cos(6*x)/24]
    @test_int [log(x/2), x, 1, -x + x*log(x/2)]
    @test_int [sqrt((1 + x)/(1 - x)), x, 3, -((1 - x)*sqrt((1 + x)/(1 - x))) + 2*atan(sqrt((1 + x)/(1 - x)))]
    @test_int [(x*log(x))/sqrt(-1 + x^2), x, 5, -sqrt(-1 + x^2) + atan(sqrt(-1 + x^2)) + sqrt(-1 + x^2)*log(x)]
    @test_int [(a + x)/(a^2 + x^2), x, 3, atan(x/a) + log(a^2 + x^2)/2]
    @test_int [sqrt(1 + x - x^2), x, 3, -((1 - 2*x)*sqrt(1 + x - x^2))/4 - (5*asin((1 - 2*x)/sqrt(5)))/8]
    @test_int [x^4/(16 + x^10), x, 2, atan(x^5/4)/20]
    @test_int [(2 + x)/(2 + x + x^2), x, 4, (3*atan((1 + 2*x)/sqrt(7)))/sqrt(7) + log(2 + x + x^2)/2]

    @test_int [x*sec(x)*tan(x), x, 2, -atanh(sin(x)) + x*sec(x)]
    @test_int [x/(-a^4 + x^4), x, 2, -atanh(x^2/a^2)/(2*a^2)]
    @test_int [1/(sqrt(x) + sqrt(1 + x)), x, 3, (-2*x^(3/2))/3 + (2*(1 + x)^(3/2))/3]
    @test_int [1/(1 - ℯ^(-x) + 2*ℯ^x), x, 4, log(1 - 2*ℯ^x)/3 - log(1 + ℯ^x)/3]
    @test_int [atan(sqrt(x))/sqrt(x), x, 2, 2*sqrt(x)*atan(sqrt(x)) - log(1 + x)]
    @test_int [log(1 + x)/x^2, x, 4, log(x) - log(1 + x) - log(1 + x)/x]
    @test_int [1/(-ℯ^x + ℯ^(3*x)), x, 3, ℯ^(-x) - atanh(ℯ^x)]
    @test_int [(1 + cos(x)^2)/(1 - cos(x)^2), x, 4, -x - 2*cot(x)]

    #= ::Section::Closed:: =#
    #=Section*7.0.?=#

    @test_int [1/(x*sqrt(-25 + 2*x)), x, 2, (2*atan(sqrt(-25 + 2*x)/5))/5]
    @test_int [sin(2*x)/sqrt(9 - cos(x)^4), x, 5, -asin(cos(x)^2/3)]
    @test_int [x^2/sqrt(5 - 4*x^2), x, 2, -(x*sqrt(5 - 4*x^2))/8 + (5*asin((2*x)/sqrt(5)))/16]
    @test_int [x^3*sin(x), x, 4, 6*x*cos(x) - x^3*cos(x) - 6*sin(x) + 3*x^2*sin(x)]
    @test_int [x*sqrt(4 + 2*x + x^2), x, 4, -((1 + x)*sqrt(4 + 2*x + x^2))/2 + (4 + 2*x + x^2)^(3/2)/3 - (3*asinh((1 + x)/sqrt(3)))/2]
    @test_int [x*(5 + x^2)^8, x, 1, (5 + x^2)^9/18]
    @test_int [cos(x)^2*sin(x)^5, x, 3, -cos(x)^3/3 + (2*cos(x)^5)/5 - cos(x)^7/7]
    @test_int [cos(4*x)/ℯ^(3*x), x, 1, (-3*cos(4*x))/(25*ℯ^(3*x)) + (4*sin(4*x))/(25*ℯ^(3*x))]
    @test_int [csc(x/2)^3, x, 2, -atanh(cos(x/2)) - cot(x/2)*csc(x/2)]
    @test_int [sqrt(-1 + 9*x^2)/x^2, x, 3, -(sqrt(-1 + 9*x^2)/x) + 3*atanh((3*x)/sqrt(-1 + 9*x^2))]
    @test_int [sqrt(4 - 3*x^2)/x, x, 4, sqrt(4 - 3*x^2) - 2*atanh(sqrt(4 - 3*x^2)/2)]
    @test_int [ℯ^(3*x)*x^2, x, 3, (2*ℯ^(3*x))/27 - (2*ℯ^(3*x)*x)/9 + (ℯ^(3*x)*x^2)/3]
    @test_int [(cos(x)*sin(x))/sqrt(1 + sin(x)), x, 3, -2*sqrt(1 + sin(x)) + (2*(1 + sin(x))^(3/2))/3]
    @test_int [x*asin(x^2), x, 3, sqrt(1 - x^4)/2 + (x^2*asin(x^2))/2]
    @test_int [x^3*asin(x^2), x, 5, (x^2*sqrt(1 - x^4))/8 - asin(x^2)/8 + (x^4*asin(x^2))/4]
    @test_int [ℯ^x*sech(ℯ^x), x, 2, atan(sinh(ℯ^x))]
    @test_int [x^2*cos(3*x), x, 3, (2*x*cos(3*x))/9 - (2*sin(3*x))/27 + (x^2*sin(3*x))/3]
    @test_int [sqrt(5 - 4*x - x^2), x, 3, ((2 + x)*sqrt(5 - 4*x - x^2))/2 - (9*asin((-2 - x)/3))/2]
    @test_int [x^5/(sqrt(2) + x^2), x, 3, -(x^2/sqrt(2)) + x^4/4 + log(sqrt(2) + x^2)]
    @test_int [sec(x)^5, x, 3, (3*atanh(sin(x)))/8 + (3*sec(x)*tan(x))/8 + (sec(x)^3*tan(x))/4]
    @test_int [sin(2*x)^6, x, 4, (5*x)/16 - (5*cos(2*x)*sin(2*x))/32 - (5*cos(2*x)*sin(2*x)^3)/48 - (cos(2*x)*sin(2*x)^5)/12]
    @test_int [cos(x)*log(sin(x))*sin(x)^2, x, 4, -sin(x)^3/9 + (log(sin(x))*sin(x)^3)/3]
    @test_int [1/(ℯ^x*(1 + 2*ℯ^x)), x, 3, -ℯ^(-x) - 2*x + 2*log(1 + 2*ℯ^x)]
    @test_int [sqrt(2 + 3*cos(x))*tan(x), x, 4, 2*sqrt(2)*atanh(sqrt(2 + 3*cos(x))/sqrt(2)) - 2*sqrt(2 + 3*cos(x))]
    @test_int [x/sqrt(-4*x + x^2), x, 3, sqrt(-4*x + x^2) + 4*atanh(x/sqrt(-4*x + x^2))]
    @test_int [cos(x)^5, x, 2, sin(x) - (2*sin(x)^3)/3 + sin(x)^5/5]
    @test_int [x^4/ℯ^x, x, 5, -24/ℯ^x - (24*x)/ℯ^x - (12*x^2)/ℯ^x - (4*x^3)/ℯ^x - x^4/ℯ^x]
    @test_int [x^4/sqrt(-2 + x^10), x, 3, atanh(x^5/sqrt(-2 + x^10))/5]
    @test_int [ℯ^x*cos(4 + 3*x), x, 1, (ℯ^x*cos(4 + 3*x))/10 + (3*ℯ^x*sin(4 + 3*x))/10]
    @test_int [ℯ^x*log(1 + ℯ^x), x, 4, -ℯ^x + (1 + ℯ^x)*log(1 + ℯ^x), -ℯ^x + log(1 + ℯ^x) + ℯ^x*log(1 + ℯ^x)]
    @test_int [x^2*atan(x), x, 4, -x^2/6 + (x^3*atan(x))/3 + log(1 + x^2)/6]
    @test_int [sqrt(-1 + ℯ^(2*x)), x, 4, sqrt(-1 + ℯ^(2*x)) - atan(sqrt(-1 + ℯ^(2*x)))]
    @test_int [ℯ^sin(x)*sin(2*x), x, 4, -2*ℯ^sin(x) + 2*ℯ^sin(x)*sin(x)]
    @test_int [x^2*sqrt(5 - x^2), x, 3, (-5*x*sqrt(5 - x^2))/8 + (x^3*sqrt(5 - x^2))/4 + (25*asin(x/sqrt(5)))/8]
    @test_int [x^2*(1 + x^3)^4, x, 1, (1 + x^3)^5/15]
    @test_int [cos(x)^3*sin(x)^3, x, 3, sin(x)^4/4 - sin(x)^6/6]
    @test_int [sec(x)^4*tan(x)^2, x, 3, tan(x)^3/3 + tan(x)^5/5]
    @test_int [x*sqrt(1 + 2*x), x, 2, -(1 + 2*x)^(3/2)/6 + (1 + 2*x)^(5/2)/10]
    @test_int [sin(x)^4, x, 3, (3*x)/8 - (3*cos(x)*sin(x))/8 - (cos(x)*sin(x)^3)/4]
    @test_int [tan(x)^3, x, 2, log(cos(x)) + tan(x)^2/2]
    @test_int [x^5*sqrt(1 + x^2), x, 3, (1 + x^2)^(3/2)/3 - (2*(1 + x^2)^(5/2))/5 + (1 + x^2)^(7/2)/7]
    end

max_stars_repo_path: test/0 Independent test suites/Stewart Problems.jl


text:
mutable struct State
    v::Int                              # Discord API version.
    session_id::String                  # Gateway session ID.
    _trace::Vector{String}              # Guilds the user is in.
    user::Union{User, Nothing}          # Bot user.
    login_presence::Dict                # Bot user's presence upon connection.
    guilds::TTL{Snowflake, AbstractGuild}     # Guild ID -> guild.
    channels::TTL{Snowflake, DiscordChannel}  # Channel ID -> channel.
    users::TTL{Snowflake, User}               # User ID -> user.
    messages::TTL{Snowflake, Message}         # Message ID -> message.
    presences::Dict{Snowflake, TTL{Snowflake, Presence}}  # Guild ID -> user ID -> presence.
    members::Dict{Snowflake, TTL{Snowflake, Member}}      # Guild ID -> member ID -> member.
    errors::Vector{Union{Dict, AbstractEvent}}            # Values which caused errors.
    lock::Threads.AbstractLock    # Internal lock.
    ttls::TTLDict                 # TTLs for creating caches without a Client.
end

State(presence::NamedTuple, ttls::TTLDict) = State(Dict(pairs(presence)), ttls)
function State(presence::Dict, ttls::TTLDict)
    presence = merge(Dict(
        "since" => nothing,
        "game" => nothing,
        "status" => PS_ONLINE,
        "afk" => false,
    ), Dict(string(k) => v for (k, v) in presence))

    return State(
        0,                          # v
        "",                         # session_id
        [],                         # _trace
        nothing,                    # user
        presence,                   # login_presence
        TTL(ttls[Guild]),           # guilds
        TTL(ttls[DiscordChannel]),  # channels
        TTL(ttls[User]),            # users
        TTL(ttls[Message]),         # messages
        Dict(),                     # presences
        Dict(),                     # members
        [],                         # errors
        Threads.SpinLock(),         # lock
        ttls,                       # ttls
    )
end

TimeToLive.TTL(s::State, ::Type{T}) where T = TTL(get(s.ttls, T, nothing))

Base.get(s::State, ::Type; kwargs...) = nothing
Base.get(s::State, ::Type{User}; kwargs...) = get(s.users, kwargs[:user], nothing)
Base.get(s::State, ::Type{Message}; kwargs...) = get(s.messages, kwargs[:message], nothing)

# Get a single channel.
function Base.get(s::State, ::Type{DiscordChannel}; kwargs...)
    return get(s.channels, kwargs[:channel], nothing)
end

# Get all the channels of a guild.
function Base.get(s::State, ::Type{Vector{DiscordChannel}}; kwargs...)
    guild = kwargs[:guild]
    haskey(s.guilds, guild) || return nothing

    return map(
        i -> s.channels[i],
        filter(i -> haskey(s.channels, i), s.guilds[guild].djl_channels),
    )
end

function Base.get(s::State, ::Type{Presence}; kwargs...)
    guild = kwargs[:guild]
    haskey(s.presences, guild) || return nothing

    return get(s.presences[guild], kwargs[:user], nothing)
end

function Base.get(s::State, ::Type{Guild}; kwargs...)
    haskey(s.guilds, kwargs[:guild]) || return nothing
    g = s.guilds[kwargs[:guild]]

    # Guilds are stored with missing channels, members, and presences.
    g = @set g.channels = map(
        ch -> s.channels[ch],
        filter(ch -> haskey(s.channels, ch), collect(g.djl_channels)),
    )
    ms = get(s.members, g.id, Dict())
    g = @set g.members = map(
        i -> get(s, Member; guild=g.id, user=i),
        filter(i -> haskey(ms, i), collect(coalesce(g.djl_users, Member[]))),
    )
    ps = get(s.presences, g.id, Dict())
    g = @set g.presences = map(
        i -> get(s, Presence; guild=g.id, user=i),
        filter(i -> haskey(ps, i), collect(coalesce(g.djl_users, Presence[]))),
    )
    g = @set g.djl_users = missing

    return g
end

function Base.get(s::State, ::Type{Member}; kwargs...)
    # Members are stored with missing user.
    haskey(s.members, kwargs[:guild]) || return nothing
    guild = s.members[kwargs[:guild]]

    haskey(guild, kwargs[:user]) || return nothing
    member = guild[kwargs[:user]]

    haskey(s.users, kwargs[:user]) || return member  # With a missing user field.
    user = s.users[kwargs[:user]]

    return @set member.user = user
end

function Base.get(s::State, ::Type{Role}; kwargs...)
    haskey(s.guilds, kwargs[:guild]) || return nothing
    roles = coalesce(s.guilds[kwargs[:guild]].roles, Role[])
    idx = findfirst(r -> r.id == kwargs[:role], roles)
    return idx === nothing ? nothing : roles[idx]
end

Base.put!(s::State, val; kwargs...) = val
Base.put!(s::State, g::UnavailableGuild; kwargs...) = insert_or_update!(s.guilds, g)
Base.put!(s::State, ms::Vector{Member}; kwargs...) = map(m -> put!(s, m; kwargs...), ms)
Base.put!(s::State, u::User; kwargs...) = insert_or_update!(s.users, u)

function Base.put!(s::State, m::Message; kwargs...)
    if ismissing(m.guild_id) && haskey(s.channels, m.channel_id)
        m = @set m.guild_id = s.channels[m.channel_id].guild_id
    end

    insert_or_update!(s.messages, m)
    touch(s.channels, m.channel_id)
    touch(s.guilds, m.guild_id)

    return m
end

function Base.put!(s::State, g::Guild; kwargs...)
    # Replace members and presences with IDs that can be looked up later.
    ms = coalesce(g.members, Member[])
    ps = coalesce(g.presences, Presence[])
    users = map(m -> m.user.id, filter(m -> m.user !== nothing, ms))
    unique!(append!(users, map(p -> p.user.id, ps)))
    g = @set g.members = missing
    g = @set g.presences = missing
    g = @set g.djl_users = Set(users)
    chs = coalesce(g.channels, DiscordChannel[])
    channels = map(ch -> ch.id, chs)
    g = @set g.djl_channels = Set(channels)
    g = @set g.channels = missing

    insert_or_update!(s.guilds, g)
    put!(s, chs; kwargs..., guild=g.id)
    foreach(m -> put!(s, m; kwargs..., guild=g.id), ms)
    foreach(p -> put!(s, p; kwargs..., guild=g.id), ps)

    return get(s, Guild; guild=g.id)
end

function Base.put!(s::State, ms::Vector{Message}; kwargs...)
    return map(m -> put!(s, m; kwargs...), ms)
end

function Base.put!(s::State, ch::DiscordChannel; kwargs...)
    if ismissing(ch.guild_id) && haskey(kwargs, :guild)
        ch = @set ch.guild_id = kwargs[:guild]
    end

    foreach(u -> put!(s, u; kwargs...), coalesce(ch.recipients, User[]))
    haskey(s.guilds, ch.guild_id) && push!(s.guilds[ch.guild_id].djl_channels, ch.id)

    return insert_or_update!(s.channels, ch)
end

function Base.put!(s::State, chs::Vector{DiscordChannel}; kwargs...)
    return map(ch -> put!(s, ch; kwargs...), chs)
end

function Base.put!(s::State, p::Presence; kwargs...)
    guild = coalesce(get(kwargs, :guild, missing), p.guild_id)
    ismissing(guild) && return p

    if haskey(s.guilds, guild)
        g = s.guilds[guild]
        g isa Guild && push!(g.djl_users, p.user.id)
    end

    haskey(s.presences, guild) || (s.presences[guild] = TTL(s, Presence))

    return insert_or_update!(s.presences[guild], p.user.id, p)
end

function Base.put!(s::State, m::Member; kwargs...)
    ismissing(m.user) && return m
    guild = kwargs[:guild]

    # Members are stored with a missing user field to save memory.
    user = m.user
    m = @set m.user = missing

    # This is a bit ugly, but basically we're updating the member as usual but then
    # filling in its user field at the same time as the user sync.
    haskey(s.members, guild) || (s.members[guild] = TTL(s, Member))
    m = insert_or_update!(s.members[guild], user.id, m)
    m = @set m.user = insert_or_update!(s.users, user)

    if haskey(s.guilds, guild)
        g = s.guilds[guild]
        g isa Guild && push!(g.djl_users, user.id)
    end

    return m
end

function Base.put!(s::State, r::Role; kwargs...)
    guild = kwargs[:guild]
    haskey(s.guilds, guild) || return r
    g = s.guilds[guild]
    g isa Guild || return r

    return if ismissing(g.roles)
        s.guilds[guild] = @set g.roles = [r]
        r
    else
        insert_or_update!(g.roles, r)
    end
end

# This handles emojis being added to a guild.
function Base.put!(s::State, es::Vector{Emoji}; kwargs...)
    guild = kwargs[:guild]

    if haskey(s.guilds, guild) && s.guilds[guild] isa Guild
        s.guilds[guild] = @set g.emojis = es
    end

    return es
end

# This handles a single emoji being added as a reaction.
function Base.put!(s::State, e::Emoji; kwargs...)
    message = kwargs[:message]
    user = kwargs[:user]
    haskey(s.messages, message) || return e

    locked(s.lock) do
        m = s.messages[message]
        isclient = !ismissing(s.user) && s.user.id == user
        if ismissing(m.reactions)
            s.messages[message] = @set m.reactions = [Reaction(1, isclient, e)]
        else
            idx = findfirst(r -> r.emoji.name == e.name, m.reactions)
            if idx === nothing
                push!(m.reactions, Reaction(1, isclient, e))
            else
                r = m.reactions[idx]
                r = @set r.count += 1
                r = @set r.me = r.me | isclient  # TODO: |= (Setfield#55).
                r = @set r.emoji = merge(r.emoji, e)
                m.reactions[idx] = r
            end
        end
    end

    return e
end

insert_or_update!(d, k, v; kwargs...) = d[k] = haskey(d, k) ? merge(d[k], v) : v
function insert_or_update!(d::Vector, k, v; key::Function=x -> x.id)
    idx = findfirst(x -> key(x) == k, d)
    return if idx === nothing
        push!(d, v)
        v
    else
        d[idx] = merge(d[idx], v)
    end
end
function insert_or_update!(d, v; key::Function=x -> x.id)
    insert_or_update!(d, key(v), v; key=key)
end


text:
# -*- coding: utf-8 -*-
# + {}
module ExpODEModels
import Base: show

using SemanticModels
using SemanticModels.ExprModels
using SemanticModels.ExprModels.Parsers
import SemanticModels: model

export ExpODEModel, model

"""    ExpODEModel

tells the model function to parse an expression as the definition
of an ODE model. Used for dispatch.
"""
struct ExpODEModel <: AbstractModel
    expr::Expr
    calls
    funcs
    variables
    domains
    values
end

""" model(::AbstractModel, expr::Expr)

dig into the expression that describes a model and break it down into components. This allows you to construct a structured representation of the modeling problem at the expression level. Just like how julia modeling frameworks build structured representations of the problems in data structures. This version builds them at the expression level.

The first argument is the type you want to construct, the second argument is the expression that you want to analyze. For example

```
model(ExpODEModel, expr)::ExpODEModel
```

"""
function model(::Type{ExpODEModel}, expr::Expr)
    matches = callsites(expr, :ODEProblem)
    @show matches
    funcs = [findfunc(expr, rhs.args[2])[1] for rhs in matches]
    params(x) = structured(x, x.args[1].args[4])
    states(x) = structured(x, x.args[1].args[3])
    fluxes(x) = structured(x, x.args[1].args[2], false)
    vars = map(x->(state=states(x),
                   flux=fluxes(x),
                   params=params(x)), funcs)
    tdomain = map(m->findassign(expr, m.args[4]), matches)
    initial = map(m->findassign(expr, m.args[3]), matches)
    return ExpODEModel(expr, matches, funcs, vars, tdomain, initial)
end

function show(io::IO, m::ExpODEModel)
    write(io, "ExpODEModel(\n  calls=$(repr(m.calls)),\n  funcs=$(repr(m.funcs)),\n  variables=$(repr(m.variables)),\n  domains=$(repr(m.domains)),  values=$(repr(m.values))\n)")
end

end


text:
using CompScienceMeshes, BEAST
using LinearAlgebra
using Profile
using LiftedMaps
using BlockArrays
using StaticArrays

ntrc = X->ntrace(X,Γ)

T = tetmeshsphere(1.0,0.3)
X = nedelecd3d(T)
Γ = boundary(T)
Y = raviartthomas(Γ)

@show numfunctions(X)
@show numfunctions(Y)

Z = BEAST.buffachristiansen2(Γ)

κ,  η  = 1.0, 1.0
κ′, η′ = √2κ, η/√2
ϵ_r =3
ϵ_b =2.0
κ_in, η_in = √6κ, η/√6

#χ = x->(1.0-1.0/ϵ_r)

function tau(x::SVector{U,T}) where {U,T}
    1.0-1.0/3.0
end

χ = tau

N = NCross()
#Volume-Volume
L,I,B = VIE.singlelayer(wavenumber=κ′, tau=χ), Identity(), VIE.boundary(wavenumber=κ′, tau=χ)
#Volume-Surface
Lt,Bt,Kt = transpose(VSIE.singlelayer(wavenumber=κ′)), transpose(VSIE.boundary(wavenumber=κ′)), transpose(VSIE.doublelayer(wavenumber=κ′))
#Kt = VSIE.doublelayerT(wavenumber=κ′)

Ls,Bs,Ks = VSIE.singlelayer(wavenumber=κ′, tau=χ), VSIE.boundary(wavenumber=κ′, tau=χ), VSIE.doublelayer(wavenumber=κ′, tau=χ)
#Surface-Surface
T  = Maxwell3D.singlelayer(wavenumber=κ)  #Outside
T′ = Maxwell3D.singlelayer(wavenumber=κ′) #Inside
K  = Maxwell3D.doublelayer(wavenumber=κ)  #Outside
K′ = Maxwell3D.doublelayer(wavenumber=κ′) #Inside

E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
H = -1/(im*κ*η)*curl(E)

e, h = (n × E) × n, (n × H) × n

@hilbertspace D j m
@hilbertspace k l o
β = 1/(ϵ_r*ϵ_b)
ν = 1/ϵ_b
α, α′ = 1/η, 1/η′
γ′ = im*η′/κ′
ζ′ = im*κ′/(ϵ_b*η′)
δ′ = im*κ′/ϵ_b

#=
eq = @varform (β*I[k,D] + η′*Lt[k,j]-γ′*Bt[ntrc(k),j] -          Kt[k,m] +
                                     (η*T+η′*T′)[l,j] -      (K+K′)[l,m] +
                                          (K+K′)[o,j] + (α*T+α′*T′)[o,m] == -e[l] - h[o])

=#

eq = @varform (β*I[k,D]-ν*L[k,D]-ν*B[ntrc(k),D] + η′*Lt[k,j]-γ′*Bt[ntrc(k),j] -          Kt[k,m] +
                    -δ′*Ls[l,D]-ν*Bs[l,ntrc(D)]            + (η*T+η′*T′)[l,j] -      (K+K′)[l,m] +
                                    -ζ′*Ks[o,D] +                 (K+K′)[o,j] + (α*T+α′*T′)[o,m] == -e[l] - h[o])

dvsie = @discretise eq  D∈X k∈X j∈Y m∈Y l∈Y o∈Y

u_n  = solve(dvsie)

#=
#preconditioner
Mxx = assemble(I,X,X)
iMxx = inv(Matrix(Mxx))

Tzz = assemble(T,Z,Z); println("dual discretisation assembled.")
Nyz = Matrix(assemble(N,Y,Z)); println("duality form assembled.")

iNyz = inv(Nyz); println("duality form inverted.")
NTN = iNyz' * Tzz * iNyz

M = zeros(Int, 3)
N = zeros(Int, 3)

M[1] = numfunctions(X)
N[1] = numfunctions(X)
M[2] = M[3] = numfunctions(Y)
N[2] = N[3] = numfunctions(Y)

U = BlockArrays.blockedrange(M)
V = BlockArrays.blockedrange(N)

precond = BEAST.ZeroMap{Float32}(U, V)

z1 = LiftedMap(iMxx,Block(1),Block(1),U,V)
z2 = LiftedMap(NTN,Block(2),Block(2),U,V)
z3 = LiftedMap(NTN,Block(3),Block(3),U,V)
precond = precond +ϵ_r*z1 + z2 + z3

A_dvsie_precond = precond*A_dvsie

#GMREs
import IterativeSolvers
cT = promote_type(eltype(A_dvsie), eltype(rhs))
x = PseudoBlockVector{cT}(undef, M)
fill!(x, 0)
x, ch = IterativeSolvers.gmres!(x, A_dvsie, Rhs, log=true,  reltol=1e-4)
fill!(x, 0)
x, ch = IterativeSolvers.gmres!(x, precond*A_dvsie, precond*Rhs, log=true,  reltol=1e-8)
=#

#Post processing
Θ, Φ = range(0.0,stop=2π,length=100), 0.0
ffpoints = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]

# Don't forget the far field comprises two contributions
ffm = potential(MWFarField3D(κ*im), ffpoints, u_n[m], Y)
ffj = potential(MWFarField3D(κ*im), ffpoints, u_n[j], Y)
ff = -η*im*κ*ffj + im*κ*cross.(ffpoints, ffm)

ffd = potential(VIE.farfield(wavenumber=κ, tau=χ), ffpoints, u_n[D], X)
ff2 = im*κ*ffd

using Plots
plot(xlabel="theta")
plot!(Θ,norm.(ff),label="far field",title="DVSIE")
plot!(Θ,√6.0*norm.(ff),label="far field",title="DVSIE 2")

import Plotly
using LinearAlgebra
fcrj, _ = facecurrents(u_n[j],Y)
fcrm, _ = facecurrents(u_n[m],Y)
vsie_j = Plotly.plot([patch(Γ, norm.(fcrj))], Plotly.Layout(title="j D-VSIE"))
vsie_m = Plotly.plot(patch(Γ, norm.(fcrm)), Plotly.Layout(title="m D-VSIE"))

#NearField
function nearfield(um,uj,Xm,Xj,κ,η,points,
    Einc=(x->point(0,0,0)),
    Hinc=(x->point(0,0,0)))

    K = BEAST.MWDoubleLayerField3D(wavenumber=κ)
    T = BEAST.MWSingleLayerField3D(wavenumber=κ)

    Em = potential(K, points, um, Xm)
    Ej = potential(T, points, uj, Xj)
    E = -Em + η*Ej + Einc.(points)

    Hm = potential(T, points, um, Xm)
    Hj = potential(K, points, uj, Xj)
    H = 1/η*Hm + Hj + Hinc.(points)

    return E, H
end

Zz = range(-1,1,length=100)
Yy = range(-1,1,length=100)
nfpoints = [point(0.0,y,z) for  z in Zz, y in Yy]

import Base.Threads: @spawn
task1 = @spawn nearfield(u_n[m],u_n[j],Y,Y,κ,η,nfpoints,E,H)
task2 = @spawn nearfield(-u_n[m],-u_n[j],Y,Y,κ_in,η_in,nfpoints)

E_ex, H_ex = fetch(task1)
E_in, H_in = fetch(task2)

Enear = BEAST.grideval(nfpoints,β.* u_n[D],X)
Enear = reshape(Enear,100,100)

contour(real.(getindex.(Enear,1)))
heatmap(Zz, Yy,  real.(getindex.(Enear,1)))

heatmap(Zz, Yy, real.(getindex.(E_in,1)))
heatmap(Zz, Yy, real.(getindex.(E_ex,1)))

contour(real.(getindex.(E_ex,1)))

Dd = range(1,100,step=1)
plot!(Yy,real.(getindex.(E_in[Dd,50],1)))
plot(collect(Yy)[2:99],real.(getindex.(Enear[2:99,50],1)))


text:
#!/usr/bin/env julia

f1 = ARGS[1]
f2 = ARGS[2]

"""
Expanding the dashed intervals into single nucleotide positions
example:
Y:2688071-2688075 =>
                    Y:2688071
                    Y:2688072
                    Y:2688073
                    Y:2688074
                    Y:2688075
"""
function expand(locci)
    list = []
    pos = split(locci, ':')
        chr = pos[1]
        region = split(pos[2], '-')
        n = parse(Int64, region[1])
        while n <= parse(Int64, region[end])
             out = "$chr" * ":" * "$n"
            push!(list, out)
            n +=  1
        end
    return list
end

"""
Reads EnsEMBL Variant Effect Predictor (VEP) Results
and returning mutation locci as a list
"""
function readVEP(fh)
    locci = []
    file  = open(fh)
    lines = readlines(file)
    for line in lines
        if ismatch(r"^#", line)
            continue
        else
            cols = split(line, '\t')
            push!(locci, cols[2])
        end
    end

    list = []
    for g in locci
        if ismatch(r"-", g)
            test = expand(g)
            for i in test
                push!(list, i)
            end
        else
            push!(list, g)
        end
    end
    return list
end

gl = readVEP(f1)
sm = readVEP(f2)

glOnly = setdiff(Set(gl), Set(sm))

"""
smOnly contains only variants excluding the common variants found in germline
"""
smOnly = setdiff(Set(sm), Set(gl))

common = intersect(Set(gl), Set(sm))
smf = open(f2)
lines = readlines(smf)
"""
fetching the variant position based annoation for the filtered
somatic variants
"""
for g in smOnly
    for line in lines
        if ismatch(r"^#", line)
            continue
        else
            cols = split(line, '\t')
            if g == cols[2]
                println(chomp(line))
            elseif ismatch(r"-", cols[2])
                pos = expand(cols[2])
                for item in pos
                    if g == item
                        println(chomp(line))
                    end
                end
            end
        end
    end
end


text:
export IPNewton

mutable struct IPNewton{T <: Real} <: AbstractOptimizationMethod{T}
    f_tol::T
    g_tol::T
    ub::Vector{T}
    lb::Vector{T}
    iterations::Int
    initial_points::Vector{Vector{T}}
end

"""
    IPNewton(model::Model, args...; kwargs...)

If an `Int` is provided the method will sample that many random initial_points and use them as
initial points for Optim IPNewton method. If a `Vector{Vector{T}}` is provided it will use them as
initial points for Optim IPNewton method.

This method provides an alternative to create box constraints. constraints can be passed as a Vector
the default constraints are `ub = Inf * ones(T, dim_unknowns(model))` and `lb = -Inf * ones(T, dim_unknowns(model))`
"""
function IPNewton(model::Model{D, T}, n_initial_points::Int; f_tol::T = T(1e-6), g_tol::T = T(1e-6),
                                                ub::Vector{T} = Inf * ones(T, dim_unknowns(model)),
                                                lb::Vector{T} = -Inf * ones(T, dim_unknowns(model)),
                                                iterations::Int = 10^5, LB::T = 0.0, UB::T = 0.6) where {D, T}

    initial_points = create_initial_points(model, n_initial_points, LB, UB)

    return IPNewton{T}(f_tol, g_tol, ub, lb, iterations, initial_points)
end

function IPNewton(model::Model{D, T}, initial_points::Vector{Vector{T}}; f_tol::T = T(1e-6), g_tol::T = T(1e-6),
                                                            ub::Vector{T} = Inf*ones(T, length(initial_points[1])),
                                                            lb::Vector{T} = -Inf*ones(T, length(initial_points[1])),
                                                            iterations::Int = 10^5) where {D, T}

    ensure_seeds_dimensions(model, initial_points)

    return IPNewton{T}(f_tol, g_tol, ub, lb, iterations, initial_points)
end

function optimize(func::Optim.TwiceDifferentiable, opt_method::IPNewton{T}, verbose::Int, i::Int) where T
    cons = Optim.TwiceDifferentiableConstraints(opt_method.lb, opt_method.ub)
    return optimize(func, opt_method, cons, Optim.IPNewton(), verbose, i)
end